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Translational  diffusion  of  water  within  biological  tissue  can  be  measured 
noninvasively  using  magnetic  resonance  imaging  techniques.  Diffusion  imaging  of 
neural  tissue  can  provide  very  valuable  information  about  the  tissue  microstructure 
because  of  its  sensitivity  to  the  molecular  motion  in  length  scales  much  smaller 
than  the  resolution  of  the  images  and  similar  to  the  dimensions  of  the  cells. 

Diffusion  attenuated  multidirectional  signal  can  be  utilized  to  provide  infor- 
mation about  the  local  orientation  of  highly  structured  environments.  However, 
the  most  common  approach  that  has  been  used  for  this  purpose  (diffusion  tensor 
imaging)  assumes  orientational  homogeneity  within  the  region.  This  assumption 
can  lead  to  incorrect  estimations  of  fiber  orientations  and  anisotropy  values.  In 
order  to  overcome  this  problem,  we  have  proposed  to  use  Cartesian  tensors  of  rank 
higher  than  2 to  model  the  measured  diffusion  coefficients.  This  approach  enables 
the  calculation  of  more  accurate  anisotropy  values  which  is  of  great  importance  in 
clinical  studies.  We  have  also  achieved  the  resolution  of  multiple  orientations  using 
higher  rank  tensors  as  well  as  through  a more  direct  approach  in  which  we  generate 
the  Laplace  series  coefficients  of  the  probability  function  on  the  surface  of  a sphere. 
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The  methods  we  present  are  applicable  to  data  that  can  be  acquired  using  clinically 
available  scanners  in  short  time  frames. 
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CHAPTER  1 
INTRODUCTION 

For  thousands  of  years,  philosophers  and  scientists  have  pondered  the  source  of 
consciousness  and  thought  and  have  believed  since  the  Middle  Ages  that  the  brain 
is  the  location  of  much  of  memory,  thought  and  emotion.  It  is  quite  surprising  that 
despite  all  the  enduring  and  widespread  interest,  the  brain,  and  the  central  nervous 
system  (CNS)  in  general,  is  the  least  understood  organ  in  human  body. 

Among  many  techniques  that  are  used  in  the  studies  involving  neural  tissue, 
the  advances  in  technology  have  made  imaging  an  indispensable  tool  in  the 
diagnosis  of  many  diseases  as  well  as  in  the  studies  aiming  a better  understanding 
of  the  brain  structure  and  function.  A remarkable  achievement  was  the  spatial 
localization  of  nuclear  magnetic  resonance  (NMR)  signal  that  eventually  led  to 
magnetic  resonance  imaging  (MRI).  Using  MRI,  it  became  possible  to  produce 
images  of  neural  tissue  with  excellent  contrast  between  different  anatomical 
structures.  Moreover,  it  was  found  that  the  NMR  signal  can  be  made  sensitive 
to  very  different  characteristics  of  the  sample,  making  it  possible  to  obtain  very 
different  information  by  using  essentially  the  same  experimental  setup.  Among 
these  characteristics  is  the  incessant  and  random  motion  of  the  molecules  being 
investigated. 

Diffusion  imaging,  by  itself,  has  many  different  applications.  Using  diffusion 
NMR,  it  is  possible  to  measure  the  diffusion  coefficients  of  many  materials. 
Furthermore,  if  the  diffusional  process  is  happening  inside  a pore,  the  dimensions 
of  the  restriction  can  be  estimated  by  using  NMR  techniques  and  the  nature  of  the 
boundaries  restricting  diffusion  can  be  understood.  In  addition,  if  the  medium  has 
an  anisotropic  structure  or  the  boundaries  are  anisotropic,  then  it  may  be  possible 
to  estimate  the  alignment  of  the  medium  or  the  orientation  of  the  boundaries. 
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1.1  The  Problem 

The  random  motion  of  water  molecules  within  the  neural  tissue  is  affected  by 
a variety  of  factors  such  as  restrictions  due  to  cell  membranes,  cytoskeleton  and 
macromolecules.  With  the  understanding  of  how  these  factors  contribute  to  the 
overall  diffusional  process,  it  may  be  possible  to  obtain  valuable  information  about 
the  biological  microstructure  simply  by  observing  the  motion  of  water  molecules. 
This  is  particularly  important  in  the  understanding  of  neural  tissue  which  is  well 
known  for  its  highly  complicated  structure. 

The  central  nervous  system  (CNS)  is  composed  of  neuronal  cells  connected  to 
each  other  with  axons  that  function  as  transmission  lines  between  different  regions. 
An  understanding  of  how  the  brain  and  the  spinal  cord  function  is  impossible 
without  the  knowledge  of  how  different  regions  are  connected  to  each  other. 
Diffusional  processes  within  the  axonal  cells  are  heavily  influenced  by  the  thin  and 
long  structure  of  these  cells.  Water  molecules  tend  to  diffuse  more  freely  along  the 
direction  of  the  fiber.  Therefore,  if  one  can  quantify  the  orientational  preference  of 
diffusion,  it  may  be  possible  to  relate  it  to  the  orientations  of  these  transmission 
lines. 

The  problem  of  mapping  the  connections  within  the  neural  tissue  resembles 
the  problem  of  creating  a road  map  of  a city  simply  by  using  the  information  about 
how  the  vehicles  are  moving  on  the  roads.  If  we  know  how  the  vehicles  are  moving 
locally  at  all  locations  we  can  create  trajectories  that  may  cover  the  whole  city. 
Moreover,  by  looking  at  certain  parameters  like  the  frequency  of  vehicles,  it  may  be 
possible  to  distinguish  the  major  highways  from  others.  Similarly  in  neural  tissue  it 
would  be  very  useful  to  have  information  about  the  integrity  and  coherence  of  the 
connections  between  several  regions  so  that  the  major  lines  can  be  differentiated 
from  others. 
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1.2  The  Method:  Diffusion  MRI 

Nuclear  magnetic  resonance  (NMR)  provides  a powerful  means  to  observe 
diffusion  noninvasively.  The  translational  motion  of  water  molecules  can  be  probed 
by  applying  magnetic  field  gradients  so  that  each  point  in  space  is  labeled  with 
a particular  magnetic  field  that  creates  a phase  change  in  the  nuclei.  A moving 
particle  experiences  different  magnetic  fields  at  different  positions.  Therefore, 
the  overall  state  of  the  molecule  is  dependent  upon  the  motional  history  of  these 
molecules.  It  has  been  observed  from  the  very  early  days  of  NMR  that  this  effect 
in  liquid  state  is  easily  detectable  by  the  NMR  spectrometer.  This  is  because  the 
timeframe  in  which  the  NMR  signal  is  produced  is  sufficient  for  the  molecules  to 
travel  a large  enough  distance  that  can  be  probed  with  relatively  small  gradients. 
Moreover,  the  resultant  signal  is  dependent  on  the  parameters  such  as  diffusion 
time  and  diffusion  gradient  strength  which  can  be  controlled  and  modified  during 
the  experiment.  Yet  another  feature  of  the  NMR  diffusion  encoding  is  that  the 
direction  along  which  these  gradients  are  applied  effectively  makes  the  signal 
sensitive  to  translational  motion  along  that  direction.  Therefore,  direction  of  these 
gradients  is  another  parameter  that  can  be  controlled. 

Repeated  measurements  of  diffusion  with  the  application  of  gradients  along 
different  directions  enable  one  to  model  the  orientational  characteristics  of  the 
diffusional  process.  If  certain  directions  are  preferred  with  respect  to  others, 
then  we  can  claim  that  those  directions  should  be  the  local  orientations  of  the 
connections.  Moreover,  if  local  orientations  from  all  locations  can  be  obtained,  that 
is,  if  the  diffusion  can  be  imaged,  then  it  may  be  possible  to  trace  the  transmission 
paths.  Finally,  if  one  can  quantify  how  much  certain  direction  is  prefered  over 
others,  he/she  may  have  information  about  the  integrity  of  the  corresponding  path. 

All  of  this  information  can  be  inferred  to  some  degree  through  diffusion 
magnetic  resonance  imaging  (MRI).  At  a coarse  level,  one  can  estimate  the  fiber 
orientations  and  quantify  the  integrity  of  these  pathways  by  measuring  anisotropy. 
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This  potentially  makes  it  possible  to  noninvasively  construct  the  anatomical,  hence 
functional,  connections  within  the  CNS.  Also,  changes  in  anisotropy  as  well  as 
orientation  of  the  axonal  cells  with  disease  provide  a potentially  useful  diagnostic 
tool. 

Despite  the  potential  of  diffusion  MRI  to  probe  the  restricted  and  anisotropic 
structures,  most  of  the  theoretical  studies  attempting  to  address  this  issue  rely 
upon  experiments  that  are  not  easy  to  implement  for  biological  samples  and 
patients.  These  requirements  are  long  acquisition  times  that  may  make  in  vivo 
imaging  impossible,  no  motion  within  the  patient  which  is  not  realistic  in  many 
situations,  and  very  strong  gradient  strengths  that  are  not  available  in  clinical 
scanners  and  may  cause  artificial  neural  stimulation  in  the  neural  tissue.  These 
factors  make  it  difficult,  if  not  impossible,  to  utilize  these  techniques  in  biological 
studies.  Therefore,  one  usually  resorts  to  phenomenological  methods  that  yield 
reliable  information  through  less  demanding  experiments. 

1.3  Outline 

In  this  dissertation,  we  present  some  of  the  developments  we  have  made  in 
the  imaging  of  diffusion  in  neural  tissue  using  magnetic  resonance  techniques. 

In  Chapter  2,  the  physics  of  diffusion  is  reviewed  briefly.  This  is  followed  by  the 
introduction  of  NMR  diffusion  experiments  that  can  be  applied  to  the  imaging 
of  biological  tissue.  Next,  we  discuss  the  first  systematic  attempt  to  quantify  the 
orientational  dependence  of  diffusional  signal  attenuation  called  diffusion  tensor 
MRI  (DT-MRI).  In  Chapter  3,  we  introduce  and  formulate  a generalization  of  this 
scheme  designed  to  overcome  some  of  the  difficulties  associated  with  DT-MRI.  This 
formalism  is  extended  to  the  formulation  of  new  rotationally  invariant  scalar  indices 
in  Chapter  4.  In  Chapter  5,  we  investigate  whether  distinct  orientations  can  be 
mapped  using  this  technique.  Also  introduced  in  this  chapter  is  another  method 
that  addresses  the  same  problem  in  a computationally  more  efficient  manner. 


CHAPTER  2 

DIFFUSION  IMAGING  FUNDAMENTALS 


2.1  Preface 

In  this  chapter  we  briefly  review  the  physics  of  diffusion  restricting  our 
emphasis  to  the  simple  case  of  free  diffusion  with  homogeneous  diffusivity.  Then 
we  discuss  the  measurement  of  diffusion  using  magnetic  resonance  techniques. 
Particularly,  we  review  the  techniques  applicable  to  the  imaging  of  the  diffusional 
properties  of  biological  tissues.  We  summarize  several  models  that  attempt  to 
explain  the  observed  signal  attenuation  with  increasing  diffusion  sensitivity. 

Finally,  we  show  how  information  inferred  using  multidirectional  diffusion  MRI  can 
be  used  to  construct  meaningful  scalar  indices  and  anatomical  connections  within 
fibrous  tissues  such  as  muscle  and  white-matter  in  the  central  nervous  system 
(CNS). 

2.2  Review  of  Diffusion  Physics 

2.2.1  Brownian  Motion 

It  was  1827  when  the  Scottish  botanist  Robert  Brown  noticed  the  perpetual 
oscillations  of  the  pollen  grains  suspended  in  water  when  he  was  working  on  the 
mechanisms  of  fertilization  in  flowering  plants  [1,  2],  This  observation  remained 
unexplained  until  Albert  Einstein,  who  was  unaware  of  Brown’s  observation  and 
searching  for  evidence  that  would  undoubtedly  imply  the  existence  of  atoms,  came 
to  the  conclusion  that  (see  Einstein  [3],  p.549)  “bodies  of  microscopically  visible 
size  suspended  in  a liquid  will  perform  movements  of  such  magnitude  that  they 
can  be  easily  observed  in  a microscope.”  1 In  his  work,  assuming  the  particles 


1 Translation  to  English  was  done  in  Fiirth  and  Cowper  [4]. 
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are  spherical,  Einstein  related  the  diffusion  coefficient  D,  to  the  fluid  viscosity  rj, 
the  gas  constant  R,  the  particles’  radius  a and  the  temperature  T through  the 
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expression 


NA  6nrja 

which  together  with  the  relation  (in  one  dimension) 


(2.1) 


(( X - X0)2)  = 2 Dt 


(2.2) 


has  enabled  the  experimental  determination  of  Avogadro’s  number,  NA  [5,  6].  The 
left  hand  side  of  the  above  expression  indicates  the  expected  value  of  the  square 
of  the  displacement  from  X0  to  X in  time  t.  Around  the  same  time  as  Einstein, 
Smoluchowski  was  able  to  reach  the  same  conclusions  using  different  means  [7]. 
2.2.2  From  Random  Walks  to  Diffusion  Equation 

A great  deal  of  theoretical  work  has  been  performed  accounting  for  the 
irregular  and  incessant  motions  exhibited  by  the  small  grains  suspended  in  a fluid. 
The  particle  undergoing  Brownian  motion  experiences  about  1021  collisions  each 
second  under  normal  conditions  [8] . Each  collision  is  thought  to  create  a sharp 
change  in  the  particle  trajectory.  Therefore  the  overall  process  has  many  degrees  of 
freedom  and  it  can  be  modeled  with  a series  of  discrete  steps  or  “random  walks.” 

If  each  step  is  identical  and  not  dependent  on  the  others,  and  the  discrete  jump 
length  has  finite  mean  and  variance,  then  as  the  number  of  steps  approach  infinity, 
the  probability  of  finding  the  particle  at  a distance  XX  away  from  its  initial 
position  after  time  Xt  is  given  by  a Gaussian  distribution, 


a consequence  of  the  central  limit  theorem  [9].  The  independence  of  the  steps 
ensure  that 


(2.3) 


’OO 


W(X  - XX,  t ) P(XX,  Xt)  d{ XX) 


(2.4) 


— OO 
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where  W ( X , t)  is  the  probability  of  finding  the  particle  at  the  position  X at  time  t. 
Expanding  the  function  W(X-  AX,  t)  around  x and  W(X,t  + At)  on  the  left  hand 
side  around  t using  Taylor  expansion  and  taking  the  limit  At  -*  0,  we  arrive  at  the 
diffusion  equation 


dW(X,t) 

~dt 


W{X,t). 


(2.5) 


Note  that  we  started  with  a particle  undergoing  Brownian  motion  under  the  in- 
fluence of  discrete  perturbations  and  reached  diffusion  equation  which  describes 
a continuous  process.  The  apparent  contradiction  is  because  diffusion  is  a macro- 
scopic manifestation  of  Brownian  motion  on  the  microscopic  level. 

For  historical  reasons  the  above  treatment  started  with  the  motion  of  small 
grains  or  any  particles  of  colloidal  size  immersed  in  a fluid.  However,  all  of  the 
reasoning  and  conclusions  presented  are  valid  if  the  observed  particle  belongs  to  the 
fluid  it  was  immersed  in.  In  this  case,  the  diffusion  process  is  called  self-diffusion. 
Throughout  this  dissertation,  we  will  be  mainly  concerned  with  the  self-diffusion 
of  water.  However,  most  of  the  time,  we  will  be  refering  to  water  self-diffusion  as 
“diffusion”  for  the  sake  of  brevity. 

2.2.3  Diffusion  in  Three  Dimensions 

In  the  above  formalism  we  have  assumed  that  the  random  walk  occurs  in  one 
dimension.  Now,  we  present  the  generalization  to  three  dimensions.  Let  us  assume 
that  we  have  a single  particle  undergoing  a series  of  random  flights  where  at  the 
J'th  step,  the  particle  moves  a distance  Ar,  = (Axj,Ayj,Azj)T  with  a probability 
distribution  function  r,(  Ax3,  Aijj,  A Zj).  Then  the  probability  of  finding  the  particle 
at  the  position 

N 

AR  = J2  Alb  (2.6) 

3=  1 

relative  to  the  initial  position  after  N steps  is  given  by  Markoff’s  method  [8]: 


/OO 

exp(— 27niq  • AR)  yfyv(q)  dq  , 

-OO 


(2.7) 


where  the  characteristic  function  is  given  by 


N 


AN{q)  ~Y\.  Ti(^rj)  exp(27rzq  • Arj)  dAr 
j= i J-°° 


j ■ 


(2.8) 


We  will  be  interested  in  the  bahaviour  of  this  general  solution  under  two 
conditions:  N — > oo  and  probability  function  Tj  is  identical  at  each  step,  i.e.,  r = Tj 
for  all  j.  If  we  further  assume  that  all  first  moments  of  displacements  vanish,  we 
get  the  expression 


^Ar(q)  = exp  ( -yqTSq  ) , 


(2.9) 


where 


' <(Ai)2> 

(AxAy) 

(AxAz) 

S : = 

(AyAx) 

((Ay)2) 

(AyAz) 

y (AzAx) 

(AzAy) 

((Az)2) 

2.2,  we  define  the  diffusion  tensor 

as 

(2.10) 


Dmn  = ±-t(ARmARn)  = ~NSn 


Substituting  this  into  Eq.  2.9  and  using  Eq.  2.7,  we  get 

1 ( ARTD_1AR\ 


(2.11) 


P(AR,  t)  = 


: exp 


V 


\J  (47rf)3det(D)  V 

This  is  the  three  dimensional  generalization  of  Eq.  2.3.  Note  that  this  expression 
is  just  an  oriented  Gaussian  function.  For  isotropic  diffusion,  the  offdiagonal 
components  of  the  diffusion  tensor  will  vanish,  and  the  diagonal  elements  will  be 
equal  to  each  other.  In  this  case,  the  probability  function  is  given  by 


(2.12) 


Piso(AR,t)  = 


1 


(4nDt)3/2 


exp  - 


|ARp 

ADt 


(2.13) 


where  D is  any  of  the  diagonal  elements  of  the  tensor.  Note  that  in  this 


case, 


(|AR|2)  = 


[ dRR 2 [ dQR 2 
Jo  Js 


(■ 4nDt )3/2 


exp 


R2 

4Dt 


= QDt 


(2.14) 
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In  the  above  expression  S is  the  unit  sphere,  R = |AR|,  and  Q is  the  solid  angle 
corresponding  to  the  direction  of  AR. 

In  three  dimensions,  by  using  a similar  analysis  as  we  have  shown  for  diffusion 
in  one  dimension,  the  differential  equation  describing  the  diffusion  process  can  be 
shown  to  be 


where  IE(R,  t)  is  the  probability  of  finding  the  particle  at  the  position  R at  time  t. 
It  is  straightforward  to  show  that  the  expression  in  Eq.  2.12  is  the  Green’s  function 
(propagator)  for  this  equation  when  no  boundary  conditions  are  specified. 

In  the  more  general  case,  when  the  probability  depends  on  the  initial  position 
of  the  particle,  we  can  define  a conditional  probability  distribution  function  that 
serves  as  the  Green’s  function;  i.e. , the  expression 


holds.  Here,  P(R|R',f)  is  the  probability  of  a particle  initially  at  the  position 
R'  to  end  up  at  R after  a time  t.  In  other  words,  P(R|R',t)  is  the  solution  to 
the  diffusion  equation,  Eq.  2.15,  when  the  particle  is  initially  at  R'.  Note  that 
as  the  system  approaches  equilibrium,  these  two  probabilities  are  expected  to  be 
essentially  the  same,  i.e., 


When  the  diffusion  process  is  restricted  by  boundaries,  then  the  solution 
depends  on  the  geometry  and  the  nature  of  these  boundaries.  The  solution  to  the 
diffusion  equation  in  a variety  of  simple  geometries  can  be  found  in  Crank  [10]. 

2.3  Measurement  of  Diffusion  Using  Magnetic  Resonance  Techniques 
Nuclear  magnetic  resonance  (NMR)  provides  a unique  way  to  quantify  the 
diffusional  characteristics  of  samples.  Because  diffusional  processes  are  dependent 
on  the  geometrical  structure  of  the  environment,  NMR  can  be  used  to  probe  the 
structural  environment  in  a noninvasive  fashion.  This  is  particularly  important  in 


(2.15) 


(2.16) 


W (R)  = lim  P(R|0,f)  . 


(2.17) 
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the  studies  involving  biological  samples  in  which  the  characteristic  length  of  the 
boundaries  restricting  diffusion  are  typically  much  smaller  than  the  resolutions 
achievable  by  conventional  magnetic  resonance  imaging  (MRI)  techniques. 

In  this  section  we  discuss  how  diffusion  affects  the  signal  detected  by  the  NMR 
spectrometer.  We  present  several  methodologies  that  are  commonly  used  to  observe 
diffusion  in  biological  samples.  We  give  emphasis  to  pulsed  field  gradient  (PFG) 
techniques  which  are  suitable  to  be  integrated  into  imaging  experiments  [11,  12]. 
Using  diffusion  MRI,  it  is  possible  to  infer  information  about  the  local  orientation, 
which  can  be  utilized  in  the  generation  of  connections  between  different  regions  of 
fibrous  tissues. 

2.3.1  Pulsed  Field  Gradient  Experiment 

A typical  nuclear  magnetic  resonance  experiment  starts  with  the  excitation 
of  the  nuclei  with  a 90°  radiofrequency  (RF)  pulse  that  tilts  the  magnetization 
vector  into  the  plane  whose  normal  is  along  the  main  magnetic  field.  Spins  that 
are  initially  coherent  dephase  due  to  many  factors,  the  most  prominent  of  which 
are  magnetic  field  inhomogeneities  and  dipolar  interactions  [13].  This  results  in  a 
decay  of  the  electromotive  force  induced  in  the  receiver.  Figure  2 1 shows  a spin 
echo  experiment  [14]  where  a subsequent  application  of  a 180°  RF  pulse  reverses 
the  dephasing  due  to  inhomogeneities  and  the  signal  is  reproduced.  The  time 
between  the  90°  pulse  and  reformation  of  the  echo  is  called  TE  and  it  is  twice 
the  time  between  the  two  RF  pulses.  This  echo  is  detected  by  a receiver  antenna 
and  is  used  to  produce  spectra.  Careful  application  of  magnetic  field  gradients 
linearly  changing  in  space  enable  the  acquisition  of  magnetic  resonance  images. 

The  details  of  the  techniques  employed  for  spatial  encoding  will  not  be  explained  in 
this  dissertation,  and  can  be  found  in  Haacke  et  al.  [15] 

An  experiment  such  as  that  depicted  in  Figure  2 1 takes  on  the  order  of  10ms 
in  a typical  imaging  experiment.  In  this  timeframe  the  spins  undergo  diffusion  that 
results  in  mixing,  which  creates  attenuation  in  the  signal  amplitude  when  magnetic 
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Figure  2 1:  Spin  echo  experiment. 

field  inhomogeneities  are  present.  This  fact  was  noticed  by  Hahn  in  the  early  days 
of  spin  echoes  [14],  but  it  was  mostly  conceived  as  an  unfortunate  phenomenon 
resulting  in  a loss  of  signal,  therefore  making  it  more  difficult  to  determine  the 
relaxation  rates.  Later  on,  Carr  and  Purcell  derived  the  echo  attenuation  when 
the  spins’  motions  are  modeled  as  a series  of  discrete  jumps  [16].  The  passage  to 
a differential  equation  was  done  by  Torrey  [17]  that  resulted  in  the  addition  of 
diffusion  terms  to  phenomenological  Bloch  equations  quantifying  the  relaxation 
phenomena.  We  will  treat  this  equation  in  detail  below. 

These  studies  assumed  that  the  magnetic  field  gradients  are  either  caused 
by  the  inhomogeneities  within  the  sample  or  an  external  magnetic  field  gradient 
is  introduced  which  is  constant  throughout  the  experiment.  In  the  1960s,  it  was 
suggested  that  this  attenuation  can  be  quantified  by  the  application  of  gradient 
pulses,  which  are  a pair  of  spatially  dependent  magnetic  fields  whose  directions  are 
along  the  main  magnetic  field.  Figure  2 2 shows  two  diffusion  gradients  added  to 
the  spin  echo  pulse  sequence  for  this  purpose.  In  this  “pulsed  gradient  spin  echo” 
(PGSE)  experiment  [18]  two  identical  gradients  around  the  180°  RF  pulse  are 
applied  with  a time  A between  their  leading  edges.  The  duration  of  these  gradients 
is  denoted  by  5.  If  G represents  the  linear  magnetic  field  gradient  applied  at  time 
ti,  after  the  application  of  the  gradient,  the  spins  at  the  position  r gain  a phase 
shift  given  by 


0 1=7 


G • r d t . 


(2.18) 
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Figure  2-2:  Pulsed  gradient  spin  echo  (PGSE)  experiment.  Two  diffusion  gradients 
G are  applied  before  and  after  the  180°  pulse. 

At  a later  time  t2,  this  phase  change  can  be  “undone”  by  the  application  of 
the  same  gradient  along  the  opposite  direction,  or  alternatively  along  the  same 
direction  but  after  the  180°  RF  pulse.  A spin  that  moves  to  a different  position 
between  the  application  of  these  pulses  experiences  a net  phase  shift  since  the 
phase  change  due  to  the  subsequent  gradient,  fc,  will  not  be  equal  to  -fa.  As  a 
result,  at  a particular  position  in  space,  there  will  be  many  particles  with  different 
phases.  The  signal  received  from  a particular  voxel  is  proportional  to  the  total 
transverse  magnetization  within  that  voxel,  given  by  the  vector  sum  of  individual 
magnetic  moments  that  will  be  denoted  by  /i: 

N 

M = ■ (2.19) 

n=  1 

ft  is  clear  that  when  the  phases  of  the  spins  (f)n  are  all  the  same,  the  magnitude 
of  the  detected  magnetization,  M,  is  equal  to  its  maximum  value  of  Nfi,  while 
when  the  phases  are  totally  random,  it  is  0.  Therefore,  the  attenuation  in  the 
magnitude  of  the  signal  indicates  the  randomness  of  the  phases  at  a particular 
position  which  in  turn  depends  on  the  randomness  in  the  spins’  motional  history 
making  it  sensitive  to  the  incoherent  motion  of  the  molecules  alone.  Diffusion 
weighted  MRI  exploits  this  phenomenon  to  quantify  diffusion  that  occurs  within 


13 


Figure  2 3:  Three  diffusion  weighted  images  from  an  excised  rat  brain  with  diffu- 
sion sensitizing  gradients  of  strength  335,  475  and  585  mT/m. 


the  sample.  Figure  2 3 depicts  the  attenuation  of  the  signal  from  an  excised  rat 
brain  with  increasing  gradient  strengths. 

Alternatively,  the  second  RF  pulse  in  a spin  echo  experiment  can  be  replaced 
by  two  90°  RF  pulses.  In  this  case,  the  generated  echo  is  called  a stimulated 
echo.  This  pulse  sequence  is  depicted  in  Figure  2-  4.  Although  there  is  a 50% 
loss  of  signal  inherent  in  this  scheme,  this  pulse  sequence  is  ideal  to  be  used  in 
experiments  where  long  diffusion  times  (A)  are  desired.  This  is  because  during  the 
time  interval  TB,  the  magnetization  is  stored  along  the  main  magnetic  field,  so  the 
spin-spin  relaxation  during  this  time  is  avoided.  Depending  on  the  relaxation  rate 
of  the  signal  and  the  length  of  TB,  when  compared  with  a corresponding  spin  echo 
experiment,  the  gain  due  to  this  may  easily  compensate  for  the  50%  signal  loss. 

I ♦ Ta  < ll * Tb  < ll »■  Ta  < 1 

90  90  90 

RF  1 

► > a < . : 

G 

— 8 

Signal  — 


Figure  2 4:  Pulsed  field  gradient  technique  using  stimulated  echoes. 
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2.3.2  From  Propagator  Formalism  to  g-Space  Imaging 

Now  we  look  at  what  happens  to  the  signal  when  the  gradient  pulses  are  much 
shorter  than  the  interval  between  them,  i.e.,  6 < A in  a PFG  experiment.  Under 
this  condition,  all  the  diffusive  processes  can  be  assumed  to  occur  in  the  absence 
of  diffusion  gradients,  and  the  diffusion  time  t can  be  taken  to  be  equal  to  A. 
Therefore,  it  is  possible  to  quantify  the  signal  attenuation  using  the  relation  [19,  20] 


where  S(G,t)  is  the  signal  detected  when  the  gradient  pulses  G are  applied  with 
a time  t in  between.  The  function  P(R|R/,t)  is,  as  before,  the  probability  for  a 
particle  initially  at  R'  to  end  up  at  R after  time  t and  p(R')  is  the  initial  spin 
density,  analogous  to  VF(R',0)  function  defined  before. 

Eq.  2.20  is  a very  important  relationship  because  it  relates  the  signal  at- 
tenuation to  the  propagator  for  the  underlying  diffusion  process.  Therefore,  if 
the  geometrical  structure  is  known,  then  one  can  calculate  the  propagator  using 
standard  techniques  and  then  relate  it  to  the  observed  signal  by  applying  Eq.  2.20. 
However,  the  inverse  problem,  i.e.,  calculating  the  propagator  from  the  observed 
signal  values,  is  in  general  not  possible  because  the  p(R')  function  is  unknown  and 
the  integration  over  R'  is  in  general  not  invertible. 

Employing  a change  of  variables  AR  = R - R'  in  Eq.  2.20,  it  is  easy  to  see 
that  the  signal  attenuations  are  given  via  an  inverse  Fourier  transform, 


/ 


dK'  p( R')  / dR  P(R|R',t)  exp  [i^5G  • (R  - R')]  , (2.20) 


/ 


(2.21) 


where 


q :=  (2tt)  , 


(2.22) 
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and  P(AR,  t)  is  the  average  propagator  [21]  defined  by 

P(AR,  t)  :=  f dR'  p(R,)P(R|R/,  t)  . (2.23) 

Employment  of  Eq.  2.21  to  calculate  the  average  propagator  is  called  “g-space” 
analysis.  Recently  it  was  shown  that  this  average  propagator  that  can  be  calcu- 
lated directly  from  the  signal  attenuations  through  a Fourier  transform  can  give 
important  information  about  compartmentation  in  biological  systems  [22]. 

Eq.  2.21  indicates  that  the  signal  attenuation  function  is  the  characteristic 
function  of  the  averaged  propagator  and  can  be  thought  of  as  the  analog  of  the 
microscopic  function  Ajv(q)  defined  in  Eqs.  2. 7-2. 8. 

2.3.3  The  Bloch- Torrey  Equation 

The  dynamics  of  the  magnetization  is  governed  by  the  phenomenological 
Bloch  Equations  [23].  Torrey  has  shown  that  when  isotropic  diffusion  is  taken  into 
account,  it  takes  the  form  [17] 

dM + 

= -ilu0M+  - ry r • GM+  - M+/T2  + P>V2M+  , (2.24) 

where  M+(:=  Mx  + iMy)  is  the  transverse  magnetization,  is  the  Larmor 
frequency,  7 is  the  gyromagnetic  ratio,  T2  is  the  spin-spin  relaxation  time  and  D 
is  the  diffusion  coefficient.  Henceforth,  we  will  call  this  equation  the  Bloch-Torrey 
equation.  2 

For  a spin  echo  experiment,  Eq.  (2.24)  can  be  solved  by  using  the  substitution 
[17] 

M+(r,  t-  G)  = m(r,  t;  G)  exp (~iu0t  - t/T2)  , (2.25) 


2 A more  complete  form  of  this  equation  will  be  treated  in  the  next  section. 
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which  puts  the  equation  in  the  form  of  the  Schrodinger  equation  for  a particle  with 
imaginary  mass  under  the  influence  of  a linear  potential  [24]: 

dm  ^ „ 

— = — zqr  • G m + DV  m . (2.26) 

For  the  free  space,  a further  simplification  is  possible  via  the  substitution 

m(r,  t\  G)  = A(t)  exp(— ryr  ■ F)  , (2.27) 

where 

F it)  :=  jT  G(0  dt-2e(t-  J 2 G(f')  dt'  , (2.28) 

and  0(x)  is  the  Heaviside  step  function  which  is  equal  to  unity  when  its  argument 
is  positive,  and  0 otherwise. 

Inserting  Eq.  2.27  into  Eq.  2.26,  we  get  a first  order  differential  equation  given 
by 

f)A 

— = -72£>|F|M  . (2.29) 

Integrating  this  equation  from  t = 0 to  t = TE  yields  the  Stejskal-Tanner  formula 
for  diffusive  attenuation  [18] 

S = S0  exp(—j252G2(A  — S/3)D) 

= S0  exp (-bD)  , (2.30) 

where  S :=  v4(TE),  S0  :=  H(0)  and 

b :=  72£2G2(A  - 5/3)  (2.31) 

is  the  6-factor.  Note  that  S is  the  magnitude  of  the  spin  echo,  and  is  the 
magnitude  of  the  signal  when  no  diffusion  gradients  are  applied. 

Enhancing  diffusive  attenuation  with  the  application  of  gradients,  as  we 
discussed  for  the  case  of  spin  echoes,  introduces  a different  contrast  mechanism 
in  magnetic  resonance  images.  The  images  of  the  signal  intensity  S , such  as 
those  in  Figure  2 3,  are  called  diffusion  weighted  images.  Diffusion-weighted 
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Figure  2-5:  The  image  with  no  diffusion  weighting  is  given  by  the  intercept  of  the 
fit  (left).  The  quantitative  diffusivity  map  is  shown  on  the  right.  Diffusion  sensitiz- 
ing gradient  is  along  the  frequency  encoding  direction  (left  to  right). 

images  have  been  utilized  extensively  in  the  imaging  of  neural  tissue  since  it  was 
shown  that  ischemic  strokes  can  be  detected  much  earlier  with  diffusion  weighted 
images  when  compared  to  traditional  Tx  and  T2  weighted  images  [25],  Moreover, 
by  repeating  the  experiment  with  different  gradient  strengths,  it  is  possible  to 
calculate  the  apparent  diffusion  coefficients  from  Eq.  2.30.  The  images  where  the 
pixels  denote  the  calculated  diffusion  coefficients  are  called  quantitative  diffusivity 
maps.  Another  quantity  that  can  be  calculated  from  Eq.  2.30  is  S0,  the  signal  if 
no  diffusion  gradient  was  applied.  Figure  2-  5 shows  the  image  with  no  diffusion 
weighting  and  the  quantitative  diffusivity  map  calculated  from  the  images  in  Figure 
2-3. 

2.4  Dependence  of  the  Signal  on  the  Gradient  Strength 
Solution  to  the  Bloch  Torrey  equation  as  demonstrated  above,  is  valid  for  free 
diffusion  as  no  boundary  conditions  were  imposed  on  the  magnetization.  However, 
diffusion  of  water  in  the  brain  is  obviously  restricted.  This  fact  creates  a deviation 
of  the  observed  signal  from  the  exponential  behavior  implied  by  Eq.  2.30.  Figure 
2 6 shows  the  magnitude  of  the  stimulated  echo  signal  intensity  observed  from 
a region  of  interest  (ROI)  specified  on  an  excised  rat  hippocampus  as  a function 
of  the  diffusion  weighting  factor  b.  The  monoexponential  behavior  is  depicted  by 
the  solid  line  in  this  logarithmic  plot.  It  is  clear  that  for  small  values  of  6,  the 
exponential  attenuation  assumption  is  valid.  In  most  applications  of  diffusion 
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Figure  2 6:  Signal  attenuation  in  arbitrary  units  as  a function  of  the  6-value  in  an 
excised  hippocampus  slice.  The  diamonds  indicate  the  experimental  observations 
of  the  signal.  The  solid  line  shows  the  exponential  attenuation  of  the  signal  im- 
plied by  the  first  three  data  points.  The  inset  shows  the  region  of  interest  (ROI) 
specified  on  the  hippocampal  slice. 

weighted  imaging  the  slope  of  this  line  is  taken  to  be  the  “apparent  diffusion 
coefficient  (ADC)”.  As  a result,  this  apparent  diffusion  coefficient  can  be  defined  to 
be 


Therefore,  the  Stejskal-Tanner  result  (Eq.  2.30)  is  applicable  to  data  with  low 
6-values,  while  as  the  diffusion  weighting  is  increased,  one  will  need  more  sophisti- 
cated models. 

Among  the  models  that  attempt  to  explain  the  nonexponential  decay  of  the 
signal,  the  biexponential  model  for  molecular  motion  in  composed  systems  [26]  has 
found  widespread  use  [27,  28].  According  to  this  model,  the  tissue  is  assumed  to  be 
composed  of  two  freely  diffusing  compartments  where  the  ADC  in  each  of  them  is 
different.  In  this  case,  the  signal  is  expected  to  decay  according  to  the  formula 


(2.32) 


S = S0  (/i  exp(-bDi)  + f2  exp(— 6D2)) 


(2.33) 
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where  Di  and  D2  are  the  ADCs  from  different  compartments,  and  /1  and  f2  are 
the  corresponding  fractions.  Note  that  since  the  signal  is  equal  to  SQ  when  6 = 0, 
the  sum  of  these  fractions  has  to  be  equal  to  1.  Therefore,  this  model  has  four 
unknown  parameters  which  can  be  taken  to  be  S0,fi,  D\  and  D2,  where 


h = 1 - /1  • 


(2.34) 


Based  on  this  model,  and  including  the  exchange  between  different  compartments, 
an  analytical  model  for  the  bovine  optic  nerve  has  been  developed  [29]. 

There  has  been  a number  of  studies  that  investigate  whether  these  diffusional 
compartments  correspond  to  intra-  and  extracellular  spaces  in  biological  samples 
[30,  31].  However,  both  theoretical  [32]  and  experimental  [31]  studies  have  shown 
that  diffusion  inside  a single  pore  can  also  give  rise  to  similar  behavior.  Another 
problem  of  the  biexponential  model  is  related  to  its  assumption  that  the  number 
of  different  compartments  is  2.  This  may  be  an  oversimplification  in  most  systems 
However,  this  model  can  easily  be  generalized  to  n such  compartments  in  which 
case  the  signal  attenuates  according  to 


S = S0  fj  exp(-bDj)  , (2.35) 

i= 1 

where 

n 

J2fi  = !•  (2'36) 

l=l 

In  order  to  find  the  value  for  n that  appropriately  describes  the  observation,  one 
needs  data  with  very  high  signal  to  noise  ratio  (SNR)  [33],  which  is  usually  not 
achievable  in  MRI  experiments  in  a reasonable  timeframe.  An  alternative  is  to  use 
the  continuous  form  of  this  model  i.e.  the  Laplace  transform:  [34] 


S = S0  / f(D)exp(—bD)dD  . (2.37) 

Jo 

Note  that  this  model  is  complicated  by  the  fact  that  the  inverse  Laplace  transform 
is  ill-conditioned. 
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Figure  2-7:  Signal  attenuation  in  arbitrary  units  as  a function  of  the  b- value  in  an 
excised  hippocampus  slice.  The  three  lines  correspond  to  three  different  fits:  Biex- 
ponential model  (solid  line),  stretched  exponential  (dotted  line),  and  Rigaut-type 
asymptotic  fractal  model  (dashed  line). 


As  an  alternative  to  multiexponential  attenuation  schemes  as  described  above, 
it  has  been  found  that  Rigaut-type  asymptotic  fractal  expression  [35,  36]  given  by 


S(b)  = S0  (l  + b/pr 


(2.38) 


well  describes  the  NMR  signal  attenuation  [37]  like  many  other  phenomena 
exhibiting  similar  concave  behavior  in  log-log  plots  [36].  Note  that  this  expression 
has  only  three  unknown  parameters:  S0,/3  and  7).  Unlike  the  multiexponential 
models  and  the  stretched  exponential  model  that  will  be  described  below,  the 
Stejskal-Tanner  type  exponential  attenuation  (as  implied  by  Eq.  2.30)  is  not  a 
special  case  of  the  expression  in  Eq.  2.38. 

Yet  another  alternative  is  to  use  the  stretched  exponential  or  Kohlrausch- 
Williams- Watts  (KWW)  function  [38,  39]: 


S{b)  = S0  exp  ( — (6DDC)Q)  , 


(2.39) 


21 


where  DDC  stands  for  “distributed  diffusion  coefficient”.  This  expression  has 
been  found  to  describe  many  physical  processes  involving  relaxation  in  disordered 
systems  when  there  is  a scale  invariant  distribution  of  relaxation  rates  [40] . 

Recently,  this  function  has  been  shown  to  yield  superior  fits  to  diffusion  data  from 
neurological  tissue  when  compared  to  the  biexponential  model  despite  having  one 
less  number  of  free  parameters  [41]. 

In  Figure  2-7  we  show  the  results  of  the  fit  of  the  data  from  excised  hippocam- 
pus to  the  biexponential,  stretched  exponential  and  Rigaut-type  asymptotic  fractal 
models.  All  models  seem  to  yield  acceptable  fits.  However,  it  is  not  possible  to 
determine  which  of  these  models  best  describes  the  diffusional  process. 

2.5  Diffusion  Tensor  Magnetic  Resonance  Imaging  (DT-MRI) 

Neural  tissue  in  central  nervous  system  (CNS)  is  primarily  composed  of  two 
distinct  areas:  white-matter  and  gray-matter  [42],  White-matter  (WM)  includes 
the  axons  that  transmit  signal  between  different  regions  of  the  CNS.  The  axons 
are  covered  with  fatty  tissue  called  the  myelin  sheath  that  isolates  the  intracellular 
space  from  extracellular  space  to  increase  the  speed  of  signal  transmission  [43]. 
Gray-matter  (GM)  is  mostly  composed  of  cell  bodies  and  nonmyelinated  axons 
where  the  axonal  architecture  is  much  more  complex. 

As  can  be  predicted  by  visually  inspecting  the  electron  microscopic  (EM) 
images  such  as  those  given  in  Figure  2 8,  the  diffusion  coefficients  along  different 
directions  in  white-matter  are  different.  The  magnetic  resonance  measurement  of 
diffusional  anisotropy  in  biological  tissue  was  first  observed  in  muscle  that  also  has 
fibrous  structure  [44] . In  white-matter  the  first  observation  of  diffusional  anisotropy 
using  MR  techniques  was  done  in  1990  [45]. 

This  diffusional  anisotropy  is  mostly  due  to  the  cellular  membranes  restricting 
the  motion  of  water  molecules,  where  myelin  and  the  cytoskeleton  are  contributing 
factors  [46].  The  black  spiral  like  structures  around  the  cell  membranes  in  the 
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Figure  2 8:  Transmission  electron  microscopic  images  from  axonal  cells.  The  left 
image  shows  a transverse  view  of  the  axons,  where  the  center  image  shows  a longi- 
tudunal  view.  The  image  on  the  right  is  a cross  section  of  the  intraaxonal  skeleton. 
Images  are  courtesy  of  Timothy  M.  Shepherd. 

transverse  view  (in  the  leftmost  image  in  Figure  2 8)  are  the  myelin  sheath.  In  the 
rightmost  image,  one  can  see  the  cytoskeleton  within  a nonmyelinated  axon. 

The  signal  attenuation  in  diffusion  weighted  MR  images  is  caused  dominantly 
by  the  diffusion  gradients  employed.  So,  the  measured  diffusion  coefficients 
depend  on  the  direction  along  which  these  gradients  are  applied.  By  repeating  the 
experiment  with  gradients  applied  along  many  directions,  it  is  possible  to  quantify 
the  diffusional  anisotropy  in  the  voxel.  This  diffusional  anisotropy  is  presumed 
to  be  related  to  the  structural  anisotropy  in  the  tissue.  The  anisotropy  maps 
produced  can  be  expected  to  give  high  contrast  between  the  highly  oriented  areas 
such  as  white-matter  and  other  regions.  Because  it  can  be  thought  of  as  a measure 
of  how  well  the  tissue  is  ordered,  an  index  quantifying  anisotropy  may  also  be 
sensitive  to  the  pathological  changes  following  many  cerebral  diseases. 

The  fact  that  diffusion  processes  are  anisotropic  can  be  used  in  yet  another 
way.  In  highly  oriented  regions,  the  water  molecules  will  preferably  diffuse  along 
the  direction  with  least  amount  of  restriction.  Figure  2 9 shows  the  effect  of 
sensitizing  the  signal  to  diffusional  processes  along  two  directions.  The  bottom 
row  shows  images  when  the  gradients  are  applied  along  the  direction  pointing 
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Figure  2-9:  Images  from  an  excised  spinal  cord  when  the  diffusion  gradients  are 
applied  with  increasing  strength  (from  left  to  right).  Top  row  includes  images  when 
the  gradients  are  perpendicular  to  the  white-matter  (WM)  fibers,  where  bottom 
row  includes  those  when  the  gradients  are  parallel  to  WM  fibers  (in  and  out  of  the 
image).  In  the  upper  leftmost  image,  GM  stands  for  gray-matter. 


perpendicular  to  the  page,  where  the  top  row  includes  the  images  when  the  gra- 
dients are  applied  along  left  and  right  directions.  When  diffusion  is  isotropic  (see 
the  free  water  surrounding  the  spinal  cord  in  Figure  2 9),  the  signal  attenuation 
seems  independent  of  the  gradient  direction.  However,  when  there  is  structural 
anisotropy,  the  change  in  the  images  with  increasing  6-values  depends  significantly 
on  the  directionality  within  the  tissue.  The  orientational  dependence  of  the  sig- 
nal attenuation  is  most  apparent  in  white-matter.  Since  water  diffusion  is  least 
restricted  along  the  fiber  directions,  the  signal  attenuates  more  rapidly  when  the 
gradients  are  applied  along  those  directions  as  seen  in  the  bottom  row.  Based  on 
this  observation,  the  direction  along  which  water  diffuses  fastest  can  be  claimed 
to  give  the  fiber  direction.  The  correspondence  between  the  fiber  orientation 
and  the  orientation  of  fastest  diffusion  is  the  main  hypothesis  behind  fiber-tract 
mapping  using  diffusion-weighted  MRI.  Starting  from  a seed  point  selected  within 
the  tissue,  repeated  stepping  in  the  direction  along  which  diffusion  is  fastest  will 
allow  the  mapping  of  fibers  passing  through  that  point  [47,  48,  49,  50].  If  the  main 
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hypothesis  of  fiber-tract  mapping  is  correct,  then  this  idea  can  be  used  to  map 
anatomically,  hence  functionally,  connected  regions  of  the  brain  and  spinal  cord. 

Estimation  of  the  Diffusion  Tensor 

Diffusion  tensor  magnetic  resonance  imaging  (DT-MRI),  or  shortly  refered  to 
as  diffusion  tensor  imaging  (DTI)  introduced  by  Basser  et  al.  [51,  52],  provides 
a simple  way  of  quantifying  anisotropy  as  well  as  estimating  the  local  fiber  direc- 
tion within  the  tissue  from  multidirectional  diffusion  MRI  data.  It  employs  the 
anisotropic  diffusion  model  previously  proposed  by  Stejskal  [20].  According  to  this 
approach,  the  diffusion  term  in  Eq.  2.24  was  replaced  by  V • (DVM+),  to  give  the 
Bloch- Torrey  equation  in  its  new  form, 

— —iwoM+  - *7 r • G M+  - M+/T2  + V • (DVM+)  - V ■ vM+  . (2.40) 

Mere  D is  a second  rank  (order),  real  valued,  positive  definite,  symmetric  tensor. 
Note  that  in  the  above  expression,  we  have  not  included  any  term  accounting  for 
the  anisotropy  of  T2.  This  is  very  well  justified  in  our  experimental  setup  because 
the  quantification  of  anisotropy  is  achieved  by  changing  the  magnitude  and  the 
orientation  of  the  diffusion  gradient  rather  than  a rotation  of  the  sample  which 
would  result  in  differing  angles  with  the  main  magnetic  field.  For  a discussion  of 
the  anisotropy  of  the  T2  relaxation  rates  in  biological  tissues,  see  Henkelman  et  al. 
[53] 

The  above  change  in  the  Bloch- Torrey  equation  from  the  form  given  in 
Eq.  2.24  yields  the  appropriate  equation  for  environments  with  organizational 
anisotropy  such  as  liquid  crystals  [54].  Here  the  last  term  is  the  convection  term 
which  was  included  for  completeness. 

Using  similar  techniques  involved  for  the  case  of  isotropic  Bloch- Torrey 
equation  given  in  Eq.  2.24,  the  full  solution  to  Bloch  equation  with  an  anisotropic 
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diffusion  term  during  the  PGSE  experiment  is  given  by 

= exP^-  luJot  -r2-nr-  [f (t)  - 20  (t  - ?£)  f (^)] 

(ItE/ 2f(t')dt'-f(W)  (*-¥)) 


D 


where  f(f)  is  defined  to  be 


f (t)  :=  / G(t')  dtf 


(2.41) 


(2.42) 


Note  that  the  last  term  in  Eq.  2.41  quantifies  the  diffusive  attenuation  and 
will  be  our  focus.  With  the  definitions  S = \M+(TE)\  and  S0  = |M+(0)|  e~TE^T21 
we  can  easily  see  that 


5^ 


exp^-72^  F(f')TDF(t')rfi'j  , 


(2.43) 


where  F (t)  was  defined  in  Eq.  2.28.  In  terms  of  f(f),  it  is  given  by 


F(t)  :=f(t)-20(f-^f  (IE 


(2.44) 


In  DT-MRI  a b-matrix,  defined  by  [51] 

r-TE 

b:=72/  F(t)F(t)Tdt, 

Jo 

is  used  to  calculate  an  effective  diffusion  tensor  through  the  equation 

5 


(2.45) 


So 


= e 


— trace(bDett) 


(2.46) 


where  it  is  assumed  that  diffusion  tensor  is  almost  constant  in  time  so  that 

Deff  ~ D. 

This  equation  has  seven  unknowns  where  one  of  them  is  the  signal  intensity  if 
there  were  no  diffusion  (50)  and  the  remaining  six  are  the  unique  elements  of  the 
diffusion  tensor.  Therefore,  seven  experiments  producing  seven  linearly  independent 
equations  are  sufficient  to  determine  the  diffusion  tensor.  Note  that  Eq.  2.46  is 
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Figure  2-10:  Diffusion  tensor  calculated  from  a series  of  diffusion 
weighted  images  of  an  excised  rat  spinal  cord.  The  order  of  the  images  is 
A cx,  Ay,  Az,  Ax,  Ay’  Az>  Ax,  Ay,  Az  from  left  to  right  and  top  to  bottom.  The 
maximum  intensity  in  the  diagonal  components  was  set  to  .00125mm2/s,  whereas 
the  offdiagonal  elements  were  scaled  between  —1.8  x 10-4  and  1.8  x 10-4mm2/s. 
The  diameter  of  the  circular  NMR  test  tube  was  5mm. 


the  corresponding  expression  to  that  in  Eq.  2.30  that  enables  the  calculation  of 
diffusion  coefficients.  Figure  2 10  shows  the  components  of  the  diffusion  tensor 
calculated  from  an  excised  rat  spinal  cord. 

We  have  detailed  two  experimental  schemes  thus  far.  In  one  of  them,  we 
can  apply  diffusion  gradients  along  a certain  direction  g with  differing  gradient 
strengths  to  calculate  the  apparent  diffusion  coefficient  D( g)  from  Eq.  2.30.  Here  g 
is  a unit  vector  representing  the  direction  of  the  diffusion  gradient,  i.e.,  G = Gg.  In 
the  second  scheme,  we  can  apply  the  diffusion  gradient  in  at  least  six  noncoplanar 
directions  and  with  differing  gradient  strengths  and  calculate  the  apparent  diffusion 
tensor  (D)  by  using  Eq.  2.46.  An  important  question  to  be  addressed  is  how  these 
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results  are  related.  If  the  background  (imaging)  gradients  are  small  compared  to 
the  diffusion  gradients,  which  is  almost  always  the  case,  then  comparison  of  Eq. 

2.30  with  Eq.  2.46  yields  the  simple  relation  [55] 

L>(g)  = gJ'Dg.  (2.47) 

Here,  we  have  implicitly  assumed  that  the  5-matrix  defined  in  Eq.  2.45  is  a pure 
outer  product  of  the  gradient  directions  scaled  by  the  5-value,  i.e., 

b ~ 5gg7  , (2.48) 

where 

5 = 47t2<3’2(A  — 5/3),  (2.49) 

as  can  be  seen  from  Eqs.  2.31  and  2.22. 

The  expression  in  Eq.  2.48  is  an  equality  (rather  than  an  approximation)  when 
the  time  dependence  of  the  gradients  is  in  their  magnitude,  i.e.,  the  orientations  of 
the  gradients  stay  the  same  throughout  the  experiment.  Therefore,  the  approxima- 
tion is  quite  accurate  when  the  diffusion  gradients  are  much  larger  than  all  other 
gradients,  which  is  true  for  most  experiments. 

Eq.  2.47  suggests  that  the  whole  diffusivity  profile  can  be  generated  from  only 
6 numbers,  namely  the  distinct  components  of  the  diffusion  tensor,  as  long  as  the 
tensor  model  well  describes  the  diffusional  process  that  occurs  within  the  sample. 
Diffusion  Tensor  Imaging  and  the  Fiber  Direction 

In  addition  to  its  great  reduction  in  the  number  of  experiments  required  to 
quantify  the  distribution  of  diffusivities  over  all  orientations,  DT-MRI  also  makes  it 
possible  to  quantify  the  prefered  direction  of  diffusion.  Note  that  Eq.  2.47  implies 
that 

S{Gg)  = S0  exp(— 5gT  D g)  . (2.50) 

We  know  that  when  the  narrow  pulse  approximation  is  valid  (i  < A),  the  average 
propagator  is  just  the  Fourier  transform  of  the  signal  attenuation  as  a result  of  Eq. 


2.21: 
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P(AR,t) 


E{ q,  t)  exp(— 27rzq  • AR)  dq  . 


Inserting  S(Gg)/S0  from  Eq.  2.50  for  E(q,t)  above,  we  get 


(2.51) 


P(AR,f) 


1 

yj  (47r£)3det(D) 


ARTD_1  AR 
4f 


(2.52) 


This  result  shows  the  consistency  of  the  Bloch- Torrey  equation  and  the  q- space 
formalism. 

The  fiber  orientations  can  be  found  by  calculating  the  directions  along  which 
this  probability  will  be  largest.  This  can  be  done  by  finding  the  peaks  of  the 
isosurfaces  of  this  function.  Since  the  only  spatial  dependence  is  in  the  exponent, 
these  isosurfaces  are  specified  by  the  relation 


AR'D^AR  = r . 


(2.53) 


Note  that  because  it  is  symmetric,  diffusion  tensor  is  diagonalizable.  Upon  diago- 
nalization,  these  isosurfaces  can  be  shown  to  be  specified  by 


m2 , on2 , (zr 

H T 1 : = T 


(2.54) 


Ai  A2  A3 

where  X',  Y1  and  Z1  are  the  components  of  AR  in  the  rotated  reference  frame, 
and  Ai,  A2  and  A3  are  the  eigenvalues  of  the  diffusion  tensor.  This  is  the  expression 
for  an  ellipsoid.  An  example  of  the  eigenvalues  and  the  corresponding  eigenvectors 
calculated  from  an  excised  rat  spinal  cord  are  shown  in  Figure  2 11. 

The  ellipsoid  generated  from  the  isosurface  of  the  probability  has  its  peak 
along  its  major  axis  defined  by  the  eigenvector  corresponding  to  the  largest  eigen- 
value i.e.,  the  principal  eigenvector.  Therefore,  the  fiber  orientation  can  be  taken 
simply  to  be  this  principal  eigenvector.  This  is  consistent  with  the  eigensystem 
shown  in  Figure  2 11  where  the  components  of  the  principal  eigenvector  are  shown 
in  the  first  row.  The  ^-component  of  the  white-matter  fibers  are  distinctively 
greater  than  the  components  in  the  image  plane. 
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Figure  2-11:  Eigenvalues  of  the  diffusion  tensor  from  an  excised  spinal  cord  (left- 
most column)  in  decreasing  order  (from  top  to  bottom).  Next  to  each  eigenvalue 
are  the  images  of  the  x,  y and  z components  of  the  corresponding  eigenvector 
(from  left  to  right).  The  maximum  intensity  in  the  eigenvalue  images  was  set  to 
0.00125mm2/s,  while  the  eigenvector  components  are  scaled  between  —1  and  1. 


Although  the  images  of  the  components  of  the  principal  eigenvector  may  be 
descriptive  in  images  like  the  spinal  cord  image  presented,  for  systems  with  more 
complicated  orientational  structure,  it  is  advantageous  to  depict  the  orientation 
in  a single  image.  This  can  be  done  by  assigning  the  components  of  the  fiber 
direction  to  the  three  different  color  channels:  red,  green  and  blue  [56].  Figure  2 12 
shows  the  orientation  maps  generated  from  the  excised  rat  brain  and  spinal  cords. 

In  these  images,  blue,  green  and  red  colors  are  assigned  to  fibers  oriented  along 
left-right,  top-bottom  directions  and  in  and  out  of  the  image  plane  respectively.  In 
these  images  the  intensity  of  the  colors  are  determined  from  the  anisotropy  of  the 
pixels. 

Scalar  Measures  Derived  From  DT-MRI 

Similar  to  what  was  presented  for  the  case  of  diffusion  weighted  imaging, 
fitting  the  data  to  the  Stejskal-Tanner  relation  for  anisotropic  diffusion,  Eq.  2.46, 
produces  an  image  with  almost  no  diffusion  weighting,  which  is  calculated  from 
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Figure  2-12:  Images  showing  the  fiber  orientations  in  excised  rat  spinal  cord  (mid- 
dle) and  brain  (right).  The  sphere  on  the  left  indicates  the  color  assignments  of 
different  orientations.  Fibers  oriented  along  the  left-right  direction  are  shown  with 
blue  colors,  where  those  fibers  oriented  along  up  and  down  direction  are  assigned 
the  green  color.  Red  color  describes  the  fibers  that  point  in  and  out  of  the  picture. 

the  intercept  of  the  fit.  The  intercept  does  not  have  any  diffusion  weighting  due 
to  the  gradients  used  in  the  pulse  sequence  since  it  corresponds  to  the  signal  at 
6 = 0.  Therefore,  the  only  diffusion  weighting  comes  from  the  magnetic  field 
inhomogeneities  within  the  sample.  This  reduction  in  the  diffusion  weighting  in  the 
MR  images  can  only  be  obtained  by  the  quantification  of  diffusion. 

A diffusion  tensor  has  9 components,  6 of  which  are  distinct  elements.  Upon 
rotations,  the  eigenvalues  will  remain  constant  where  the  components  of  the 
eigenvectors  will  change.  Therefore,  an  (rotationally)  invariant  basis  of  the  diffusion 
tensor  has  three  elements  implying  that  one  can  write  three  independent  indices 
to  describe  the  rotationally  invariant  features  of  the  geometry  associated  with  the 
tensor.  These  three  elements  can  be  chosen  to  be  the  three  eigenvalues.  However,  it 
is  not  easy  to  interpret  the  meaning  of  the  eigenvalues  perhaps  except  the  principal 
(greatest)  eigenvalue,  which  is  just  the  diffusivity  along  the  fiber  direction.  Another 
choice  could  be  the  first  three  moments  of  the  diffusion  tensor,  namely  trace(D), 
trace(D2)  and  trace(D3).  All  of  these  indices,  like  the  eigenvalues,  have  units.  In 
the  context  of  DT-MRI,  the  first  moment  of  the  diffusion  tensor,  giving  the  average 
of  diffusivities  along  all  directions,  has  been  chosen  to  be  the  first  scalar  measure 
derived  from  the  diffusion  tensor.  This  quantity,  called  mean  diffusivity  [57]  is 
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given  by 


(D)  = 


trace(D)  A^  + A2  -h  A3 


(2.55) 


3 


3 


Traditionally,  it  has  been  found  useful  to  have  unitless  scalars  instead  of 
the  second  and  third  moments.  Especially  an  index  that  measures  the  level 
of  organization  in  a pixel  would  be  very  useful  since  it  would  make  it  easy  to 
distinguish  highly  organized  regions,  like  white-matter,  from  others.  There  have 
been  several  suggestions  for  the  parametrization  of  anisotropy  that  may  accomplish 
this  goal  in  the  tissue.  The  most  commonly  used  anisotropy  index  is  Fractional 
Anisotropy  (FA)  [57],  which  is  expressed  in  terms  of  the  eigenvalues  as 


1 in  the  complete  anisotropic  case.  As  for  the  third  invariant,  a scalar  measure 
such  as  skewness  [58]  can  be  used.  However,  in  neural  tissue  the  third  invariants 
have  not  found  widespread  utilization.  In  Figure  2 13  we  show  the  nondiffusion 
weighted  image  of  the  rat  brain  calculated  from  the  intercept  of  the  fit  along  with 
mean  diffusivity  and  fractional  anisotropy  maps. 

Fiber-Tract  Mapping  in  Neural  Tissue  using  DT-MRI 

The  orientation  information  obtained  from  the  diffusion  tensor  as  described 
above  can  be  utilized  in  the  construction  of  maps  of  connections  between  different 
regions  of  the  fibrous  tissues.  This  may  be  accomplished  by  repeatedly  stepping 
in  the  direction  along  the  prefered  direction  of  diffusion,  which  is  given  by  the 
principal  eigenvector.  The  constructed  fiber-tracts  can  be  thought  of  as  curves 
whose  tangent  vector  is  set  equal  to  this  direction.  The  equation  of  motion  in  this 
case  is  given  by  a Frenet  formula  [49]: 


(2.56) 


FA  takes  the  value  of  0 for  totally  isotropic  tensors  whereas  it  takes  the  value  of 


dr 


(2.57) 
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Figure  2 13:  Nondiffusion  weighted  image  S0  (left  column),  mean  diffusivity  (mid- 
dle column),  and  fractional  anisotropy  (right  column)  images  from  the  excised  rat 
spinal  cord  (top  row)  and  brain  (bottom  row). 


where  r is  the  position  vector  within  the  image,  ds  is  the  infinitesimal  scalar 
distance  on  the  curve,  and  e>  is  the  principal  eigenvector  from  diffusion  data. 

This  equation  can  be  solved  by  integration  with  user  supplied  initial  condition. 
Therefore  the  initial  condition  specifies  the  points  where  the  curves  originate. 

An  important  issue  in  this  process  is  to  decide  when  the  tract  tracing  should 
be  terminated.  As  the  tracts  enter  regions  with  isotropic  structure,  the  fiber 
direction  becomes  uncertain.  Therefore,  an  anisotropy  index  can  serve  as  a 
termination  criterion  in  fiber-tract  mapping.  Figure  2-14  shows  an  example  of 
calculated  fibers  from  normal  and  injured  spinal  cords  when  the  ROIs  are  selected 
on  two  planes  in  the  dorsal  column,  and  the  calculation  was  done  by  integrating 
Eq.  2.57  using  fourth  order  Runge-Kutta  technique  [59].  Tracking  is  terminated 
when  fractional  anisotropy  drops  below  a prespecified  value  of  0.25.  As  a result  the 
injury  site  suffers  a discontinuity  in  the  fibers. 

Figure  2 15  shows  two  major  white-matter  pathways  calculated  in  an  excised 
rat  brain.  The  seed  points  were  selected  in  corpus  callosum  and  cerebral  peduncles 
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Figure  2-14:  White  matter  fiber  tracts  of  normal  (left)  and  injured  (right)  rat 
spinal  cords  are  shown  in  blue.  The  seed  points  were  selected  in  the  dorsal  columns 
of  the  respective  spinal  cords.  The  tract  images  were  overlaid  on  images  with  no 
diffusion  weighting. 


and  fibers  were  tracked  in  orthograde  as  well  as  retrograde  directions.  Tract  tracing 
is  terminated  when  the  FA  value  drops  below  0.175. 

Problems  of  DT-MRI  Based  Fiber  Tracking 

Despite  the  fact  that  only  limited  quantitative  verification  of  the  DT-MRI 
based  fiber-tract  mapping  results  has  been  done,  qualitative  comparisons  with 
known  anatomy  of  the  nervous  tissue  has  shown  that  DT-MRI  is  able  to  correctly 


Figure  2-15:  White  matter  fiber  tracts  calculated  in  the  cerebral  peduncles  (left) 
and  corpus  callosum  (right)  in  an  excised  rat  brain.  The  fiber  tracts  (in  blue)  are 
overlaid  on  fractional  anisotropy  maps.  These  brain  samples  were  placed  in  15 mm 
NMR  tubes. 
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map  the  major  axonal  pathways  [50].  Despite  the  promising  results  achieved  by 
DT-MRI,  it  is  known  that  it  has  significant  weaknesses. 

One  of  the  problems  with  DT-MRI  is  that  the  formalism  presented  to  calculate 
the  diffusion  tensor  from  signal  attenuation  does  not  account  for  spatial  dependence 
of  the  diffusion  tensor  within  a voxel.  It  is  the  correct  description  for  a medium 
that  is  homogeneously  anisotropic  such  as  liquid  crystals.  It  is  not  known  how 
spatial  dependence  of  the  tensor  would  affect  the  measured  signal  attenuation. 

Another  problem  with  the  formalism  just  presented,  is  that  it  assumes  no 
boundary  conditions,  hence  assumes  free  diffusion  and  does  not  account  for 
restricting  boundaries.  These  problems  are  present  in  diffusion  weighted  imaging 
as  well.  Therefore,  the  word  ‘apparent’  has  been  used  before  the  phrases  ‘diffusion 
coefficient’  and  ‘diffusion  tensor’  to  emphasize  that  calculated  diffusivities  are  in 
reality  dependent  on  the  parameters  of  the  pulse  sequence  used. 

DT-MRI  based  fiber  tracking  assumes  that  there  is  a single  fiber  direction 
within  a voxel.  A typical  axon  has  a diameter  of  about  5 fim  and  typical  voxel 
volume  of  the  images  we  acquire  at  our  institution  is  on  the  order  of  (100 pm)3.  So, 
in  every  voxel,  there  is  a bundle  of  axons.  It  is  known  that  in  many  areas  of  the 
brain,  these  axons  cross.  Fiber  crossing  within  a voxel  may  result  in  deviations  in 
the  calculated  directions  and  premature  termination  of  tracking  due  to  a decrease 
in  the  anisotropy  value.  This  effect  will  be  more  serious  in  clinical  imaging,  where 
voxel  volume  will  be  in  the  order  of  1 mm3.  Similar  problems  may  occur  even 
when  there  is  a single  fiber  bundle  in  the  voxel  with  a curvature  on  the  order 
of  1/d  or  higher,  where  d is  the  length  of  the  edge  of  a voxel.  Overcoming  the 
difficulties  associated  with  voxels  with  orientational  heterogeneity  will  be  the  topic 
of  upcoming  chapters. 


CHAPTER  3 

GENERALIZATION  OF  DIFFUSION  TENSOR  MR, I 
3.1  Introduction 

Until  recently,  quantification  of  diffusional  characteristics  in  anisotropic  bio- 
logical tissues  has  mainly  employed  the  diffusion  tensor  MRI  formalism  described 
in  the  previous  chapter.  This  model  was  originally  developed  to  measure  the  dif- 
fusional motion  of  molecules  that  have  organizational  anisotropy,  such  as  in  liquid 
crystals  [20].  However  unlike  liquid  crystals,  anisotropy  in  tissue  is  mainly  due  to 
geometric  restrictions  to  water  translational  motion.  Diffusion  processes  inside 
the  tissue  are  extremely  complex  making  accurate  mathematical  modeling  very 
difficult  [27,  29,  60],  and  experiments  that  test  these  models  are  very  demanding 
[61,  30].  Even  though  diffusion  tensor  imaging  applications  have  used  a relatively 
simple  rank-2  tensor  model  of  diffusion,  very  informative  maps  of  anisotropy  as  well 
as  fiber  directions  have  been  generated  that  make  fiber  tract  mapping  possible  in 
highly  structured  tissues  [47,  48,  49,  50]. 

Despite  the  early  success  of  diffusion  imaging,  it  has  been  shown  [62]  that  the 
rank-2  diffusion  tensor  model  has  important  limitations;  the  most  important  of 
which  occurs  when  there  is  orientational  heterogeneity  of  the  fibers  in  the  voxel 
of  interest.  To  overcome  this  difficulty,  Tuch  et  al.  [63]  have  proposed  the  use  of 
diffusion  imaging  with  diffusion-weighting  gradients  applied  along  many  directions 
distributed  almost  isotropically  on  the  surface  of  a unit  sphere.  This  “high  angular 
resolution  diffusion  imaging  (HARDI)”  method  does  not  assume  any  a priori 
knowledge  on  the  diffusivity  profile  like  that  assumed  in  the  rank-2  diffusion 
tensor  model.  In  this  scheme,  apparent  diffusion  coefficient  along  each  direction  is 
calculated  using  the  isotropic  version  of  the  Stejskal-Tanner  equation,  Eq.  2.30. 
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Recently,  Frank  introduced  an  improvement  in  the  HARDI  method  [64]  that 
transforms  the  measured  distribution  of  diffusivities  into  components  of  a spherical 
tensor  of  high  ranks.1  In  this  method,  the  spherical  harmonic  transform  (SHT)  of 
the  diffusivity  profile  is  calculated  using 


is  truncated  so  that  only  the  most  significant  terms  are  included  in  the  expansion. 
However,  it  is  not  clear  how  this  scheme  can  be  integrated  into  the  Bloch- Torrey 
formalism.  For  example,  the  series  truncation  indicates  that  the  diffusion  coeffi- 
cients, as  calculated  from  the  isotropic  Bloch- Torrey  equation,  Eq.  2.24,  must  be 
altered.  However  this  approach  is  not  consistent  with  the  Bloch- Torrey  equation 
with  Stejskal’s  modification  for  anisotropic  media,  Eq.  2.40,  because  the  truncated 
series  might  include  terms  higher  than  rank-2.  Therefore  a more  general  description 
is  needed  to  relate  the  measured  diffusion  profile  to  the  model  used  to  describe  the 
diffusion  process. 

In  this  chapter,  we  propose  an  extension  to  the  formal  description  of  spin 
diffusion  to  include  Cartesian  representation  of  higher  rank  (larger  than  2)  dif- 
fusion tensors.  The  Bloch-Torrey  equation  [17]  is  restated  and  solved  to  yield  an 
expression  corresponding  to  the  Stejskal- Tanner  relation  [18]  that  quantifies  the 
diffusive  signal  attenuation.  This  approach  helps  resolve  the  ambiguity  in  the 
process  of  calculating  the  spherical  harmonic  transform  of  the  diffusivity  profile  and 
the  termination  of  the  Laplace  series.  Moreover  using  the  Stejskal-Tanner  relation, 


1 The  words  “rank”  and  “order”  are  used  interchangeably  in  tensor  analysis  liter- 
ature. In  the  present  work  we  will  use  the  former. 


(3.1) 


The  resulting  Laplace  series 


OO 


(3.2) 


37 


generalized  to  higher  rank  tensors,  allows  the  straightforward  calculation  of  all  the 
components  of  higher  rank  diffusion  tensors  by  using  a least  squares  fitting  routine. 
This  makes  it  unnecessary  to  evaluate  the  spherical  harmonic  transform  from  the 
diffusivity  profiles,  which  is  a computationally  difficult  task. 

We  show  how  the  relations  between  the  coefficients  of  the  Laplace  series  and 
components  of  the  generalized  diffusion  tensors  can  be  calculated  analytically. 

The  termination  of  the  Laplace  series  expansion  at  a certain  order  lmax  (effectively 
setting  all  higher  rank  components  to  zero)  is  analogous  to  fitting  the  data  to  a 
rank-/max  Cartesian  tensor  model.  Once  a Cartesian  tensor  of  a particular  rank 
is  calculated,  the  components  of  a lower  rank  tensor  can  be  determined  by  using 
simple  linear  relations  making  it  unnecessary  to  refit  the  data  to  a lower  rank 
tensor  model.  For  convenience,  these  relationships  between  tensor  components  up 
to  rank-6  are  provided.  Also  we  derive  the  relations  between  the  diffusivities  and  a 
rank-2  diffusion  tensor,  even  when  the  diffusion  profile  has  higher  rank  components. 
We  show  by  simulation  the  inadequacy  of  using  the  traditional  rank-2  tensor 
model  and  show  the  diffusion  profiles  using  a higher-rank  tensor  model  of  various 
anatomical  regions  on  an  excised  rat  brain  image  calculated  using  generalized 
diffusion  tensor  imaging  methods. 

3.2  Generalized  Diffusion  Tensor  Imaging 

We  have  shown  in  the  previous  chapter  that  when  the  narrow  pulse  approxi- 
mation is  met,  the  observed  signal  attenuation  is  just  the  characteristic  function  of 
a macroscopic  (average)  displacement  probability  function  (see  Eq.  2.21).  In  the 
microscopic  domain,  in  three  dimensions,  the  characteristic  function  obeyed  Eq. 

2.9.  It  was  shown  in  Eq.  2.50  that  when  the  direction  of  the  gradients  had  no  time 
dependence,  i.e.,  when  the  condition 


G(f)*gG(f) 


(3.3) 
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is  met,  DT-MRI  assumes  a macroscopic  characteristic  function  that  has  essentially 


the  same  form  with  its  microscopic  counterpart.  As  a result,  the  displacement 
probability  functions  as  implied  by  DT-MRI  are  given  by  Gaussian  functions  whose 
isosurfaces  are  ellipsoids.  However,  experiments  have  suggested  that,  although  it 
is  possible  to  calculate  the  local  directionality  when  there  is  a single  orientation 
within  the  volume,  when  there  is  orientational  heterogeneity,  the  ellipsoids  of 
DT-MRI  yield  incorrect  results. 

In  order  to  model  more  complicated  angular  structures,  we  have  proposed  a 
generalization  to  DT-MRI  in  which  the  characteristic  function  of  the  underlying 
displacement  probability  has  been  assumed  to  be  given  by 

E(Gg)  = exp  I -b  ■ • ■ ^2  Diih...i,9ii9i2  ■ • ■ ) , (3.4) 


Z2=l 


U= l 


where  Dili2 are  the  components  of  the  Cartesian,  rank-/  tensor,  and  g is  a 
unit  vector  whose  components  can  be  written  in  terms  of  the  polar  angle  9 and 
azimuthal  angle  </>  as 
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^ cos  9 i 

(3.5) 


In  order  to  achieve  the  signal  attenuation  as  described  in  Eq.  3.4,  we  start 
with  an  extension  to  the  Bloch- Torrey  equation  [17]  to  include  a phenomenological 
diffusion  term  with  a high-rank  (>  rank  2)  Cartesian  tensor: 


dM+ 

dt 


—iio0M+  — i'y r ■ GM+  — M+/T2 

3 3 3 

+EE-E  9i\9i 2 ■ ■ • 9ii^  M+  . 

*1=1  *2=1  *1=1 


(3.6) 


Upon  the  substitutions  of  Eqs.  2.25  and  2.27  into  Eq.  3.6,  one  gets  a gen- 
eralized Stejskal-Tanner  formula,  when  the  background  (imaging)  gradients  are 
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negligible  compared  to  the  diffusion  gradients,  given  by 
In  S = In  So  — b x 

3 3 3 

'y  ] y v • ■ y ' ^hh-.-kdiiSh  . . . gtl  . (3-7) 

il=l  i2=l  ii=l 

This  equation  makes  it  possible  to  calculate  all  the  components  of  the  diffusion 
tensor  of  general  rank  by  means  of  a simple  multilinear  regression.  Comparison 
of  Eq.  2.30  and  Eq.  3.7  implies  that  the  diffusion  coefficient  along  the  gradient 
direction,  as  included  in  Eq.  2.24  will  be  given  by 

3 3 3 

77(g)  = EE  -E  Diii2. ,.ii9ii9i2  • ■ • 9ii  • (3.8) 

i i=lJ2  = l ii  = 1 

If  l is  an  odd  number,  then  this  relation  implies  that 

77(-g)  = -77(g),  l is  odd.  (3.9) 

However,  since  negative  diffusion  coefficients  are  nonphysical,  l is  forced  to  be  an 
even  number.  Therefore  from  Eq.  3.8,  we  arrive  at  the  condition  for  antipodal 
symmetry  of  the  diffusivities, 

D{- g)  = 77(g),  l is  even  . (3.10) 

A general  rank-/  Cartesian  tensor  has  3;  terms,  which  is  a very  large  number 
for  higher  ranks.  For  example,  a rank- 10  tensor  will  have  59049  components. 
However,  symmetries  provide  a very  significant  reduction  in  the  number  of  distinct 
components.  This  follows  from  the  realization  that  77^^  ^ is  a totally  symmetric 
tensor.  Total  symmetry  is  due  to  the  fact  that  this  tensor  links  the  components  of 
the  same  vector  to  a scalar  (77(g)).  For  example,  in  the  case  of  l = 2, 

3 3 3 3 

77(g)  = 77  ij  g3gi 

i=  1 j= 1 i~l  j= 1 

3 3 3 3 

— ^]l  9*9]  — 'y~'l  'y^,  77 jj  gigj 

j= 1 i=l  i=l  j= 1 


(3.11) 
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implies  Dtj  — Dji  since  it  is  true  for  all  vectors  g.  A similar  analysis  for  the  general 
/ case  yields 


where  (*ii2  ■ ■ - ii)  stands  for  all  permutations  of  the  indices.  This  symmetry  reduces 
the  number  of  distinct  elements  to  [65] 


which  is  only  66  for  / = 10  case. 

To  use  Eq.  3.7  to  derive  the  distinct  components  of  the  rank-/  diffusion  tensor, 
we  will  need  to  know  how  many  times  a given  element  is  repeated.  We  will  call 
this  the  multiplicity  of  that  element,  and  denote  it  with  the  letter  fj,.  Knowing  the 
multiplicity  of  every  unique  element,  we  can  rewrite  Eq.  3.7  as 


where  Dk  is  the  k-th  unique  element  of  the  tensor,  and  gk(p)  is  the  component  of 
the  gradient  direction  specified  by  the  p- th  index  of  k-th  unique  element  of  the 
generalized  diffusion  tensor.  The  multiplicity  of  a component  of  a rank-/  tensor  is 
given  by 


where  nx,  ny  and  nz  are  respectively  the  number  of  x , y and  z indices  included 
in  the  full  sequence  of  subscripts  defining  the  component  of  the  tensor.  As  an 
illustration,  the  unique  components  of  a generalized  diffusion  tensor  up  to  rank  8 
and  their  multiplicities  are  listed  in  Table  3-1. 

3.2.1  Relationships  Between  Laplace  Series  of  HARDI  and  Traditional  (/  = 2)  DTI 
In  this  section  it  is  assumed  that  a second  rank  tensor  is  sufficient  to  correctly 
characterize  the  diffusion  coefficients  along  all  directions,  i.e. , Eq.  2.47  is  correct, 


(3.12) 


(/  + 1)(Z  4-  2) 


2 


(3.13) 


(3.14) 


(3.15) 


Table  3 1:  Distinct  components  of  the  diffusion  tensor  up  to  rank-8,  and  their 
multiplicities. 
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to  find  the  corresponding  coefficients  in  the  Laplace  series.  Then  Eq.  2.47  is 
substituted  into  Eq.  3.1  and  the  integrals  evaluated  to  give  the  relationships 


0-00  = 
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The  same  procedure  can  be  repeated  to  calculate  the  correspondence  between  the 
coefficients  in  the  spherical  harmonics  expansion  and  Cartesian  tensors  of  higher 
rank  by  using  Eq.  3.8  and  Eq.  3.1. 

3.2.2  Effect  of  Higher  Rank  Components  on  the  Construction  of  Diffusion  Tensor 
Now,  we  will  assume  that  the  distribution  of  diffusivities  is  correctly  charac- 
terized by  a tensor  of  arbitrary  rank  and  investigate  what  happens  when  we  fit  the 
data  to  a rank-2  tensor  model.  To  do  this,  we  assume  that  data  is  acquired  using  a 
very  efficient  scheme,  as  in  Frank  [64] , where  a single  measurement  is  performed  at 
6 = 0,  and  many  measurements  along  different  directions  (g)  are  performed  with  a 
finite  (and  constant)  b-value.  In  this  case,  diffusivity  along  g is  given  by 

D(  g)  = -Jln^.  (3.17) 

0 Do 

If  we  fit  the  same  data  to  a rank-2  tensor  model  by  using  multilinear  regression, 
then  components  of  the  diffusion  tensor  will  be  such  that 

/ dg  (ln  f + 6gT]Dg)  (3.18) 

is  minimized.  The  minimum  value  is  achieved  when  (using  Eq.  3.17) 


[ cfe(gTDg)2  -2  f dgL>(g)grDg 

a J a 


= o . 


(3.19) 


Here  5 denotes  the  variations  of  the  quantity  in  brackets  with  small  changes  in 
the  tensor  components.  If  we  define  the  first  term  in  square  brackets  as  Ix  and  the 
second  term  as  —2 /2,  then  performing  the  relevant  simple  integrations  yields 
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(3.22) 


We  can  solve  Eq.  3.19  for  the  elements  of  the  diffusion  tensor  by  setting  the 
partial  derivatives  equal  to  zero: 

d(h  - 2/2) 
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(3.23) 


For  the  diagonal  elements,  this  yields  the  matrix  equation 
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where  for  the  offdiagonal  elements  of  the  diffusion  tensor  we  have 


Dij  — /jj  , for  i ^ j . 


(3.25) 


These  are  the  expressions  relating  a general  distribution  of  diffusivities  to  the 
calculated  components  of  a rank-2  tensor.  Note  that  in  the  derivations  we  did  not 
assume  that  diffusivities  are  characterized  by  pure  quadratic  forms  of  a rank-2 
tensor.  Therefore,  Eqs.  3.24  and  3.25  are  the  generalized  versions  of  Eq.  2.47. 

Now  we  insert  the  spherical  harmonic  decomposition  of  the  diffusivities  into 
Eq.  3.22  and  calculate  the  components  of  the  diffusion  tensor  from  Eqs.  3.24 
and  3.25.  To  evaluate  the  integrals,  we  expand  g^gj  in  a Laplace  series  and  find 
that  only  components  up  to  / = 2 are  nonzero.  Then,  orthogonality  of  spherical 
harmonics  makes  it  simple  to  find  the  components  of  the  diffusion  tensor  in  terms 


of  the  components  of  the  spherical  tensor: 
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(3.26) 


A close  examination  of  these  equations  show  that  these  are  exactly  equivalent  to 
Eqs.  3.16,  which  were  derived  under  the  assumption  that  diffusivities  along  all 
directions  are  pure  quadratic  forms,  i.e.,  with  no  / > 2 terms. 

The  equivalence  of  Eqs.  3.26  with  Eqs.  3.16  proves  that  even  if  there  are 
higher  rank  components  in  the  distribution  of  diffusivities,  these  components  do 
not  corrupt  the  results  achieved  by  fitting  the  data  to  a rank-2  tensor  model.  In 
other  words,  if  the  calculations  are  done  by  evaluating  the  components  of  the 
Laplace  series  and  setting  the  components  with  / > 2 equal  to  zero,  then  the  results 
are  equivalent  to  those  achieved  by  fitting  the  data  to  a Cartesian  rank-2  tensor 
model.  Furthermore,  since  the  equivalence  of  the  two  approaches  stems  from  the 
orthogonality  of  the  spherical  harmonics  as  outlined  above,  it  is  also  true  for  fitting 
the  distribution  of  diffusivities  to  a general  rank-/  Cartesian  tensor  model. 

3.2.3  Components  of  a Lower  Rank  Tensor  from  Those  of  a Higher  Rank  Tensor 

Once  the  data  is  fitted  to  a rank-/  tensor  model,  it  is  possible  to  calculate  the 
components  of  the  lower  rank  tensors  without  refitting  the  data  to  the  lower-rank 
tensor  model.  Calculation  of  the  components  of  lower  rank  tensors  from  those  of 
higher  rank  tensors  is  done  by  employing  analytical  linear  relationships  between  the 
components  of  lower  and  higher  rank  tensors.  Here,  we  show  the  derivation  of  the 
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Figure  3 1:  Diffusion  ellipsoids  and  diffusivity  profiles  with  increasing  anisotropy. 

relation  between  the  0-th  rank  tensor  to  the  components  of  the  second  rank  tensor. 
Then  we  will  give  the  results  of  the  similar  derivations  for  tensors  of  rank  up  to  6. 

We  start  with  the  relation  between  the  rank-0  tensor  and  the  corresponding 
component  of  the  spherical  tensor  that  is  calculated  by  using  Eq.  3.1.  In  the  case 
of  rank-0,  there  is  only  one  component  of  the  Cartesian  tensor  which  will  be  labeled 
as  D.  Using  Eq.  3.1,  it  is  easy  to  see  that 

D = aoo/v^  , (3.27) 

and  all  other  components  will  vanish. 

Next,  we  look  at  the  expressions  relating  the  higher  rank  Cartesian  tensor 
components  to  the  components  of  the  Laplace  series.  These  expressions  are  already 
given  in  Eqs.  3.16  for  a rank-2  tensor.  Inserting  the  expression  for  a0 0 into  Eq. 

3.27,  we  get  the  desired  relation 

D = ^(A<x  + Dyy  + Dzz)  . (3.28) 

As  expected,  this  expression  is  the  same  as  that  for  ‘mean  diffusivity’  (see  Eq. 

2.55).  In  a similar  fashion,  relations  between  tensors  of  other  ranks  can  be  derived. 
We  provide  these  relations  up  to  tensors  of  rank  6 in  Table  3-2. 
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Table  3 2:  Expressions  relating  the  components  of  Cartesian  tensors  up  to  rank  6 
to  those  of  lower  rank  tensors. 


D — |(DxX  + Dyy  + D 


Dxx  — 35(9  D 
^vv  = 35  $D 

Dzz  = 3 

-Dxv  = f ( Dxxxv  + D 


XXXX  T 8DXXyy  + 8D 

X 


yyyy  “h  8Dxxyy  T 8D. 


yyzz 


yy  — 35 

35(9DZZZZ  + 8DXXZZ  + 8Dyyzz 

xy  ^(.•‘-^xxxy  i J-^yyyx  i ^zzxy  ) 


Dyyyy  Dzzzz  2Dyyzz) 

CXxxx  Dzzzz  2Z)xxzz) 


Dx 


Dyyyy  %Dxxy  y) 


+ -DZZ) 


•Dxz  — f(Dxxxz  + A;zzx  + Dyyxz) 

DyZ  — -j(Dyyyz  + Dzzzy  + Dxxyz) 


wv 


Dyyyy  — 


Dzzzz 


Dxxxy  — 
Dxxxz 
Dyyyx  = 
Dyyyz  = 
Dzzzx 
Dzzzy  = 
Dxxyy  — 


Dyyzz 


Dxxyz 
Dyyxz  — 
D^zxy 


231  (^3-DxXXXXX  + Dyyyyyy  + D ZZZZZZ  + 2 4 Z)xxxxyy  + 2 4 -DXXXXZZ  + 3 DyyyyZZ 

+3  DZZZZyy  — 18Dyyyyxx  — 18Z)zzzzxx  — 36  -Dxxyyzz) 

23l(43.Dyyyyyy  + -Dxxxxxx  + Dzzzzzz  + 24Dyyyyxx  + 24  -Dyyyyzz  + 3 -Dxxxxzz 
+3-Dzzzzxx  18Dxxxxyy  - 18Dzzzzyy  - 36 Dxxyyzz) 

231  (43-DZZZZZZ  + T>XXXXXX  + Dyyyyyy  + 24  D zzzzxx  + 24DZZZZyy  + 3 D XXXXyy 
T3  Dyyyyxx  — l8Dxxxxzz  — 1 8 Dy  y y y zz  — 36  DXXyyZZ) 

22  ( 3)DXXXXXy  + 4j9xxxyyy  + 4 DXXXZZy  — Dyyyyyy,  — D ZZZZXy  ~ 2DyyyZZX) 

22  ( ^^xxxxxz  + 4DXXXZZZ  + 4Dxxxyyz  — Dzzzzzx  — DyyyyXZ  — 2 Dzzzyyx ) 

22  ( ^Dyyyyyx  + 4/JXxxyyy  + 4Z2yyyzzx  D XXXXXy  D ZZZZXy  2D  XXXZZy^ 

22  ( 5DyyyyyZ  + 4 Dy  y y zzz  + 4Dyyyxxz  ~ DZZZZZy  ~ Dxxxxyz  — 2D  ZZZXXy) 

22  ( ^DzZZZZX  T 4 DXXXZZZ  + 4 DZZZyyX  DXXXXXZ  DyyyyXZ  2Z9XXXyyz) 

22  ( ^ DzZZZZy  + '1  DyyyZZZ  + 4 DZZZXXy  — DyyyyyZ  — Dxxxxyz  — 2Dy 

1386  ( 321  DXXXXyy  + 321  Dyyyyxx  + 306DXXyyZZ  — 36DXXXXZZ  — 36tx'yyyyZZ 

1 9Z^XXXXXX  ^9  Dyy  yyyy  15D 

zzzzxx  15Dzzzzyy  + 2 D 

ZZZZZZ ) 

1386  (321Dxxxxzz  + 321DZZZZXX  + 306Dxxyyzz  — 36  Dxxxxyy  — 36Dzzzzyy 

19DXXXXXX  19DZZZZZZ  15DyyyyXX  15Dyyyyzz  + 2 Dyyyyyy  ) 

1386  (321  Dyyyyzz  + 32 1 D zzzzyy  + 306Dxxyyzz  — 36Dyyyyxx  — 36  Dzzzzxx 

— 19Dyyyyyy  — 19DZZZZZZ  — 1 5 DXXXXZZ  — 15DXXXXyy  + 2 D XXXXXX) 

66  ( 1 17 dxxxxyz  + 16Dyyyxxz  + 16Dzzzxxy  — DyyyyyZ  — Dzzzzzy  — 2 Dyyyzzz) 
66  ( 1 7 Dyyyyxz  + 16DXXXyyz  + 16DZZZyyX  — DXXXXXZ  — D ZZZZZX  — 2 D XXXZZZ) 

( 1 f ^ZZZZXy  T l6DXXXZZy  + 16DyyyZZX  DXXXXXy  ~ ^Vyyyyx  2 / ^xx  xyyy  ) 


j ) 

-'yyyxxzj 

■ 36Dy 


66 


47 


3.3  Results 

Traditionally,  visualization  of  the  rank-2  diffusion  tensors  has  been  done  by 
constructing  the  ellipsoid  given  in  Eq.  2.54.  Note  that  the  axes  of  the  resulting 
surfaces  will  have  the  units  of  length.  However,  higher  rank  tensors  do  not  have 
matrix  representations.  Therefore,  it  is  not  easy  to  generalize  linear  algebraic 
operations,  such  as  diagonalization  and  inversion,  to  higher  rank  tensors.  Instead  in 
the  visualization  of  HARDI  data,  a parametrized  2-surface  given  by  [66] 

— {D(9,4>)  sin  9 cos  0,  D(9,  </>)  sin  9 sin  cf),D(9,(f))  cos  9)  (3.29) 

has  been  used  to  depict  the  diffusivity  profile.  In  this  case,  the  units  of  the  axes 
are  that  of  diffusivity.  As  an  illustration  of  the  correspondence  between  these  two 
surfaces,  Figure  3-  1 shows  the  calculated  surfaces  when  diffusivities  are  described 
by  rank-2  diffusion  tensors  of  increasing  anisotropy  where  the  principal  axis  is  not 
changed.  It  is  obvious  that  the  same  directional  information  is  contained  in  the 
two  schemes,  where  increased  anisotropy  is  characterized  by  thinning  of  the  surface 
near  origin  in  the  latter  case.  Visualization  of  the  diffusion  tensors  by  using  the 
parametrized  surface  given  by  Eq.  3.29  makes  it  possible  to  easily  generalize  the 
visualization  to  diffusivities  specified  by  higher  rank  tensors.  This  is  accomplished 
by  inserting  Eq.  3.8  into  Eq.  3.29. 

When  a general  diffusivity  profile  is  fitted  to  a rank-2  tensor  model,  a signifi- 
cant part  of  the  information  contained  in  the  diffusivity  profile  can  be  lost.  This  is 
evident  from  Table  3-2,  which  shows  that  the  components  of  a lower  rank  tensor  is 
just  a weighted  average  of  the  components  of  the  higher  rank  tensor.  To  illustrate 
this  loss  of  information,  four  rank-4  diffusion  tensors  were  simulated  so  that  they 
all  correspond  to  the  same  rank-2  tensor  through  the  relationships  in  Table  3-2.  On 
the  left  side  of  Figure  3 -2  are  these  four  possibilities  of  diffusivity  profiles  produced 
by  choosing  four  different  rank-4  tensors.  The  parametrized  surface  corresponding 
to  the  rank-2  tensor  that  corresponds  to  all  of  these  rank-4  tensors  is  given  on  the 
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Figure  3 2:  On  the  left  are  four  diffusion  profiles  from  four  different  rank-4  tensors 
chosen  so  that  they  will  all  correspond  to  the  same  rank-2  tensor  through  the  rela- 
tionships given  in  Table  3-2.  The  right  figure  shows  the  diffusivity  profile  for  this 
rank-2  tensor. 


right.  Note  that  only  in  the  first  case  on  the  left  are  the  diffusivities  implied  by  the 
rank-2  tensor  model  consistent  with  the  actual  rank-4  tensor.  The  other  three  plots 
clearly  show  the  inadequacy  of  using  a lower  rank  tensor  model. 

To  provide  an  experimental  test  of  the  high-rank  Cartesian  tensor  modeling, 
we  acquired  a high  angular  resolution  diffusion  weighted  image  of  an  excised  rat 
brain.  The  components  of  a Cartesian  tensor  of  a particular  rank  were  calculated 
using  Eq.  3.14  using  linear  regression.  In  Figure  3 3,  we  show  the  unique  elements 
of  the  diffusion  tensors  up  to  rank  6 for  the  same  slice  of  the  excised  rat  brain 
image.  These  tensor  component  images  are  ordered  in  the  same  order  as  the 
components  in  Table  3-1.  In  order  to  provide  a relatively  uniform  presentation, 
the  maximum  values  in  the  images  are  scaled  to  4.0  x 10-4,  5.0  x 10-4,  6.0  x 10~4 
and  7.0  x 10 ~4mm2/s  for  ranks  6,  4,  2 and  0 respectively.  It  is  apparent  that  all 
the  components  that  have  even  number  of  x’s,  y’s  and  z’s  in  the  subscripts  are 
distinctively  greater  than  the  other  components.  For  selected  regions  of  interest 
in  the  brain,  we  have  calculated  the  diffusivity  profiles  corresponding  to  Cartesian 
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RANK=6 


Figure  3 3:  Distinct  components  of  Cartesian  diffusion  tensors  of  ranks  6,  4,  2 and 
0 for  an  excised  rat  brain.  The  diameter  of  the  circular  NMR  test  tube  used  is 
15mm. 
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tensors  of  increasing  rank  up  to  rank  8.  The  results  are  shown  in  Figure  3 4 with 
an  image  in  the  first  column  representing  the  fractional  anisotropy  [57]  from  the 
rank-2  tensor.  The  diffusivity  profiles  of  increasing  rank  are  shown  in  the  second 
through  sixth  column  for  each  region  of  interest  defined  in  the  fractional  anisotropy 
image.  The  slice  positions  in  Rows  a,  d,  and  e are  the  same  and  correspond 
approximately  to  position  of  the  slice  in  Fig.  40  in  Paxinos  and  Watson  [67].  Also 
the  slice  position  in  Rows  b,  c and  f correspond  approximately  to  position  of  the 
slices  in  Figs.  60,  49  and  38  respectively  of  Paxinos  and  Watson  [67].  The  regions 
of  interest  in  each  slice  are  indicated  with  a circle  with  an  attached  handle  in  each 
fractional  anisotropy  image. 

3.4  Discussion 

The  formalism  we  have  introduced  is  not  a theoretical  solution  for  the  case 
of  restricted  diffusion  that  occurs  in  biological  tissues,  but  rather  a phenomeno- 
logical model  that  attempts  to  quantify  the  diffusivities  calculated  from  HARDI 
experiments  in  terms  of  a Cartesian  tensor  of  any  rank.  From  a conceptual  point 
of  view,  it  is  equivalent  to  constructing  the  Laplace  series  of  the  spherical  harmonic 
transform  for  the  diffusivities.  Therefore,  the  modified  Bloch- Torrey  equation  can 
be  seen  as  the  appropriate  generalization  for  highly  structured  materials  and  ac- 
tually the  equation  that  is  implied  when  the  diffusion  profile  is  represented  by  the 
results  of  a spherical  harmonic  transform.  The  equivalence  of  the  results  obtained 
from  the  generalized  tensor  approach  and  spherical  harmonic  transform  makes  it 
unnecessary  to  calculate  the  coefficients  in  the  Laplace  series  since  the  Cartesian 
tensor  possesses  the  same  information  and  can  be  calculated  using  simpler  regres- 
sion methods.  Moreover,  since  all  components  of  the  generalized  diffusion  tensor 
have  real  values,  post-processing  schemes  employing  the  Cartesian  representation  of 
the  tensor  are  simpler  to  implement. 

When  a diffusion  coefficient  measurement  is  performed  in  the  presence  of 
restrictions,  by  employing  the  Stejskal-Tanner  method,  the  calculated  diffusion 
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Figure  3-4:  Diffusivity  profiles  from  different  anatomical  regions  and  their  evolu- 
tion as  the  rank  of  the  model  is  increased.  Each  row  is  taken  from  a voxel  within 
the  following  anatomical  region:  a)  Corpus  callosum,  b)  Facial  nerve,  c)  Lateral  en- 
torhinal  cortex,  d)  Medial  lemniscus,  e)  Dentate  gyrus  and  f)  Fasciculus  retroflexus. 
These  ROIs  are  shown  with  a circle  with  an  attached  handle  on  FA  maps  in  the 
leftmost  column. 


coefficient  represents  the  apparent  diffusion  rate  and  in  reality  is  dependent  on 
the  parameters  of  the  pulse  sequence.  This  inconsistency  of  observed  diffusion 
coefficients  from  experiments  performed  with  different  parameters  is  a consequence 
of  the  fact  that  the  model  assumes  free  and  homogeneous  diffusion  within  the 
voxel  and  hence  it  is  severely  under- parameterized.  Traditional  rank- 2 diffusion 
tensor  imaging  makes  the  same  assumptions  but  also  assumes  that  the  diffusing 
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medium  has  anisotropy,  which  was  originally  used  to  describe  the  diffusion  of  very 
long  and  thin  molecules  when  their  orientations  are  not  random  [54],  In  most 
applications  of  diffusion  tensor  imaging,  the  magnetic  resonance  signal  comes  from 
water  molecules  that  clearly  do  not  fit  this  description.  In  an  abstract  fashion, 
the  ellipsoidal  diffusion  of  water  molecules  in  fibrous  tissue  can  be  viewed  as 
representing  a process  similar  to  the  diffusion  of  very  long  and  thin  molecules  in  a 
molecular  environment  like  liquid  crystals.  In  this  situation,  the  diffusion  profile 
in  a voxel  of  tissue  may  be  thought  of  as  representing  diffusion  of  water  in  a single 
pipe  in  such  a way  that  the  water  molecules  behave  as  a whole  like  the  molecules 
in  a liquid  crystal,  i.e.  the  water  diffuses  mainly  along  the  pipe.  This  situation  may 
be  visualized  in  a physically  meaningful  way  with  a rank-2  tensor  as  an  ellipsoidal 
displacement  profile  or  a “peanut”  shaped  diffusivity  profile,  as  shown  in  part  (a)  of 
Fig.  3 4. 

However  water  diffusion  in  tissue  may  have  much  more  complex  structure,  such 
as  might  be  represented  by  multiple  intersecting  “pipes”  along  which  the  water 
diffuses.  This  situation  may  not  be  adequately  modeled  with  rank-2  tensors  in 
which  case  the  higher-rank  tensor  model  becomes  more  appropriate.  In  such  com- 
plex situations,  the  higher-rank  tensor  model  of  the  observed  diffusion  profile  does 
not  have  a direct  physical  analogy.  Rather,  this  model  provides  a mathematical 
description  of  the  diffusion  profile  that  allows  visualization  using  the  parametrized 
2-surface  given  in  Eq.  3.29.  Although  using  a scheme  that  does  not  take  restricted 
diffusion  into  account  may  be  a problem  from  a theoretical  perspective,  it  is  the 
strength  of  this  approach  from  an  experimental  point  of  view,  since  it  is  possible  to 
create  contrast  using  relatively  small  gradients  and  in  short  time  frames.  Even  with 
these  limitations,  the  value  of  information  gained  makes  this  method  a useful  tool 
for  the  quantification  of  anisotropy  and  potentially  in  the  determination  of  fiber  di- 
rections. For  example  since  the  appearance  of  the  profile  from  the  corpus  callosum 
in  part  (a)  of  Fig.  3 4 does  not  change  as  the  rank  of  the  model  increases  above 
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rank  2,  the  profiles  may  represent  a voxel  with  all  the  fibers  oriented  in  the  same 
direction.  As  shown  in  Figure  3-4,  the  profiles  of  part  (b)  and  (f)  are  similar  in  ap- 
pearance at  all  ranks,  except  for  an  overall  rotation,  and  imply  more  heterogeneity 
in  tissue  structure  than  is  implied  by  profiles  in  part  (a).  The  same  can  be  stated 
for  the  profiles  in  part  (c),  (d)  and  (e).  Although  fine  structure  differences  in  the 
profiles  are  visible,  it  is  difficult  to  make  predictions  with  this  limited  data  set.  In 
addition,  the  profile  differences  visible  here  might  be  enhanced  if  the  measurements 
were  extended  to  higher  6-values  and  the  overall  fine  structure  in  the  diffusivity 
profiles  above  rank  4 or  6 shown  here  may  be  the  result  of  noise  in  the  data. 

The  practical  limit  to  the  rank  of  the  tensor  model  is  given  by  the  number  of 
directions  along  which  the  data  is  collected,  since  the  number  of  unique  compo- 
nents of  the  Cartesian  tensor  given  by  Eq.  3.13  must  be  less  than  or  equal  to  the 
number  of  directions  along  which  the  diffusivity  is  measured.  Also  the  rank  of  the 
tensor  determines  the  number  of  distinct  fiber  directions  that  can  be  resolved  using 
a particular  acquisition  scheme.  However,  this  is  not  a serious  limitation  because 
one  would  rarely  expect  to  resolve  more  than  three  or  four  different  orientations 
within  a voxel.  Therefore  any  deviation  in  diffusivity  profiles  observed  with  higher 
rank  tensor  components  is  likely  to  be  due  to  noise,  so  our  approach  provides  a 
way  to  reduce  the  noise  in  the  fitted  diffusivity  profiles  for  lower  rank  components. 
Note  that  as  the  rank  of  the  tensor  model  is  increased,  the  data  becomes  more 
prone  to  effects  of  noise  since  the  number  of  parameters  deduced  from  the  same 
data  is  increased.  In  addition,  it  is  possible  to  define  certain  measures  to  decide 
which  tensor  model  best  describes  the  diffusivity  within  a voxel,  as  was  done  for 
the  spherical  harmonic  transform  approach  [68],  Then  using  the  relations  between 
the  components  of  higher  and  lower  rank  tensors,  a voxel  can  be  represented  with  a 
dynamically  selected  tensor  model. 

Although  diffusivity  profiles  are  able  to  indicate  the  complexity  of  the  fiber 
structure  within  the  voxel,  it  is  not  clear  how  one  can  extract  the  distinct  fiber 
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orientations  from  them.  It  was  shown  that  the  maxima  of  the  diffusivity  profiles 
may  not  correspond  to  distinct  fiber  directions  [69].  A possible  alternative  would 
be  to  make  a similar  analysis  for  the  level  surfaces  of  the  averaged  propagator  as 
implied  by  generalized  diffusion  tensor  imaging  and  HARDI.  We  will  be  addressing 
this  issue  in  Chapter  5. 

3.5  Conclusions 

We  have  presented  an  extension  to  the  Bloch-Torrey  formalism  for  the 
calculation  of  diffusion  profiles  in  tissue  that  employs  Cartesian  tensors  of  rank 
higher  than  the  usual  rank-2  tensor  model.  The  corresponding  Stejskal-Tanner 
formula  allowed  the  direct  calculation  of  the  components  of  the  higher  rank  tensor. 
Our  method  yields  the  same  results  as  the  previously  proposed  spherical  harmonic 
expansion  of  the  diffusivities  without  any  need  to  evaluate  the  spherical  harmonic 
transform.  We  have  shown  the  inadequacy  of  fitting  the  data  to  a rank-2  tensor 
model  both  by  simulations  and  on  images  from  an  excised  rat  brain  with  modest 
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CHAPTER  4 

GENERALIZED  SCALAR  INDICES 
4.1  Preface 


This  chapter  details  the  derivation  of  rotationally  invariant  scalar  measures 
from  higher  rank  diffusion  tensors,  as  well  as  from  functions  defined  on  unit  sphere. 
This  is  accomplished  by  using  an  expression  that  generalizes  the  evaluation  of  the 
trace  operator  to  tensors  of  arbitrary  rank  and  even  to  functions  whose  domains 
are  the  unit  sphere.  We  show  that  mean  diffusivity  is  invariant  to  the  selection 
of  tensor  rank  for  the  model  being  used.  However,  this  rank  invariance  does  not 
apply  to  the  anisotropy  measures.  Therefore,  we  propose  a variance  based,  general 
anisotropy  measure.  Also  an  information  theoretical  parametrization  of  anisotropy 
is  introduced  that  is  frequently  more  consistent  with  the  meaning  attributed  to 
anisotropy.  We  accomplish  this  by  associating  anisotropy  with  the  amount  of 
orientational  information  present  in  the  data  regardless  of  the  imaging  technique. 
Using  a simplified  model  of  fibrous  tissue,  we  simulate  anisotropy  values  with 
varying  orientational  complexity  and  tensor  models.  Simulations  suggest  that  a 
lower  rank  tensor  model  may  produce  artificially  low  anisotropy  values  in  voxels 
with  complex  structure.  This  is  confirmed  with  a spin  echo  experiment  performed 
on  an  excised  rat  brain. 

4.2  Introduction 

Diffusional  attenuation  of  the  magnetic  resonance  (MR)  signal  as  a result  of 
the  mixing  of  phase  incoherent  spins  has  been  known  since  Hahn  introduced  spin 
echoes  in  the  early  days  of  MR  [14].  We  have  shown  that  this  attenuation  due 
to  diffusion  can  be  enhanced  by  the  application  of  so-called  diffusion  gradients 
enabling  the  measurement  of  diffusion  constants.  For  samples  with  anisotropic 
structure,  this  attenuation  depends  on  the  orientation  along  which  these  gradients 
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are  applied.  Therefore,  by  repeating  the  diffusion  measurement  with  diffusion 
gradients  applied  along  different  directions,  it  is  possible  to  measure  the  diffusional 
anisotropy  which  is  presumed  to  be  related  to  structural  anisotropy  within  highly 
structured  environments,  such  as  muscle  and  white-matter  in  neural  tissue. 
Quantification  of  anisotropy  provides  a previously  inaccessible  contrast  mechanism 
in  MRI.  It  is  of  great  significance  to  investigate  whether  this  contrast  mechanism 
may  provide  previously  unavailable  diagnosis  of  any  pathological  conditions  or 
whether  it  can  complement  the  existing  diagnostic  techniques.  For  a review  of  the 
sensitivity  of  these  scalar  indices  to  different  pathologies,  see  Dong  et  al.  [70]. 

Many  of  the  previously  introduced  scalar  indices  assume  that  the  model  being 
used  is  traditional  (rank-2)  DTI  [52,  71,  72,  73,  74,  58].  The  failure  of  traditional 
DTI  in  the  presence  of  orientational  heterogeneities  may  create  problems  not 
only  in  the  determination  of  fiber  orientations,  but  also  in  the  calculated  scalar 
measures.  Figure  4 1 shows  simulations  of  a region  with  crossing  fibers.  It  is 
evident  from  the  second  image  that  an  anisotropy  index  based  on  rank-2  DTI,  such 
as  fractional  anisotropy  (FA)  [57],  produces  significantly  low  values  when  there 
are  fiber  crossings.  This  is  partly  due  to  the  fact  that  the  orientational  variation 
in  the  diffusivities  is  less  in  these  regions.  However,  traditional  DTI  suffers  from 
a further  reduction  of  anisotropy  values  as  a result  of  the  excessive  smoothing 
introduced  by  the  employment  of  the  rank-2  tensor.  This  is  apparent  in  the 
diffusivity  profiles  implied  by  the  rank-2  and  rank-6  diffusion  tensors  as  presented 
in  the  two  rightmost  images  of  Fig.  4 1.  Since  most  of  the  clinical  studies  quantify 
DTI  through  mean  diffusivity  and  a measure  of  anisotropy,  a reexamination  of  the 
derivation  of  these  indices  is  of  great  importance. 

Previously  the  variance  of  diffusivities,  measured  along  different  directions 
with  the  HARDI  method,  was  proposed  as  an  anisotropy  index  [75]  that  doesn’t 
assume  the  rank-2  tensor  model.  This  approach  has  a number  of  problems.  First 
of  all  the  calculated  “anisotropy”  maps  have  the  same  units  as  diffusivity  and 
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Figure  4-1:  Simulations  of  an  environment  with  crossing  fibers.  Left  image  shows 
the  fiber  bundles  simulated.  The  resulting  FA  image  implied  by  the  rank-2  tensor 
model  is  shown  in  the  second  figure.  The  last  two  figures  show  the  diffusivity  pro- 
files from  selected  pixels,  in  the  image  matrix  of  256  x 256,  implied  by  rank-2  and 
rank-6  tensors  respectively. 

the  images  produced  are  diffusivity  weighted.  This  may  create  problems  in  the 
interpretation  of  contrast  (and  lack  of  contrast)  seen  in  the  images  because  the 
observed  variations  may  be  due  to  variations  of  diffusivity  and  not  just  anisotropy. 
Second,  the  range  of  values  this  index  can  take  is  unclear,  thus  making  it  difficult 
to  scale  the  images  in  a consistent  way.  Also,  since  this  approach  uses  only  the 
discrete  samples  of  the  diffusivity  profile,  the  computed  values  depend  on  the 
distributions  of  the  gradient  vectors  on  the  unit  sphere.  (Note  that  this  distribution 
is  never  truly  isotropic  except  when  the  directions  are  specified  by  Platonic  solids 
[76].)  Finally  since  this  measure  of  anisotropy  is  derived  from  discrete  samples  with 
no  functional  fit,  one  may  expect  this  index  to  be  very  sensitive  to  noise. 
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In  this  chapter,  we  revisit  the  problem  of  quantification  of  mean  diffusivity 
and  anisotropy  with  the  purpose  of  formulating  measures  generalized  to  higher 
rank  tensors  as  well  as  to  functions  whose  domains  are  the  unit  sphere.  We 
show  that  the  commonly  used  expression  for  mean  diffusivity  (D)  is  a model 
independent  measure  and  simply  equal  to  the  rank-0  tensor.  This  is  not  true  for 
the  anisotropy  measures  such  as  FA  and  relative  anisotropy  (RA)  that  are  functions 
of  the  variance  of  the  eigenvalues  of  the  diffusion  tensor,  which  is  not  equal  to 
the  variance  of  the  diffusivities  along  all  directions.  Therefore,  we  construct  a 
generalized  anisotropy  (GA)  index  that  is  based  on  the  variance  of  the  normalized 
diffusion  coefficients.  The  normalization  step  utilizes  a generalized  expression 
for  the  trace  operation  and  removes  the  undesired  diffusion  weighting  from  the 
resulting  images.  Similarly,  we  propose  another  index  based  on  the  entropy  of  the 
function  defined  on  the  unit  sphere.  This  information  theoretical  formulation  of 
anisotropy,  which  we  called  the  scaled  entropy  (SE),  treats  the  normalized  function 
as  a probability  distribution  function.  The  construction  of  both  SE  and  GA  indices 
ensure  that  the  resulting  values  are  in  the  interval  [0, 1),  0 corresponding  to  the 
isotropic  profile.  We  provide  exact  expressions  for  (D),  and  GA  indices  for  tensors 
up  to  rank-6,  and  present  images  of  these  measures  for  a HARDI  data  set  from  an 
excised  rat  brain. 

Using  simulations  of  a simplified  model  of  fibrous  tissue,  we  show  that 
anisotropy  measures  calculated  from  a rank-2  tensor  model  may  be  significantly 
smaller  than  the  anisotropy  values  calculated  from  higher  rank  tensors,  when  there 
is  more  than  one  fiber  direction.  We  show  several  regions  in  our  data  in  which  this 
effect  is  significant. 

We  discuss  whether  a simple  calculation  of  an  information  theoretical 
anisotropy  index,  for  the  case  of  a rank-2  tensor,  is  possible  and  formulate  an 
index  called  visual  anisotropy  (VA)  that  is  defined  in  terms  of  the  von  Neumann 
entropy.  Just  like  RA  and  FA  can  be  thought  of  as  the  simple  implementation  of 
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a variance  based  anisotropy  index  in  the  case  of  a rank-2  tensor,  VA  is  the  cor- 
responding simple  implementation  of  an  entropy  based  index  for  rank-2.  Finally, 
we  present  several  information  theoretical  scalar  measures  that  may  be  useful  in 
the  comparison  of  two  functions  whose  domains  are  the  unit  sphere.  Although 
we  present  our  results  using  diffusivity  profiles,  the  formulations  remain  valid  for 
any  other  positive  valued  data  acquired  on  or  projected  on  to  the  surface  of  a unit 
sphere. 


4.3.1  Generalization  of  the  Trace  Operator 

The  trace  of  a rank-2  tensor  is  given  by  the  sum  of  diagonal  elements  of 
the  tensor  when  it  is  represented  by  a matrix.  It  is  straightforward  to  show  that 
trace  is  rotationally  invariant.  Trace  can  also  be  expressed  as  the  integral  of 
the  quadratic  forms  of  the  tensor  which  enables  the  generalization  of  the  trace 
operation  to  higher-rank  tensors  and  even  to  functions  defined  on  the  surface  of  a 
sphere: 


Here,  g is  a unit  vector,  A is  any  rank-2  tensor  and  S is  the  unit  sphere.  Note  that 
since  the  integrand  in  the  above  expression  has  antipodal  symmetry,  Eq.  4.1  is  also 
valid  when  the  integral  is  evaluated  on  the  unit  hemisphere  that  will  be  denoted  by 
Q.  In  this  case,  the  coefficient  before  the  integration  should  be  replaced  by  3/27 r: 


It  is  possible  to  generalize  this  operation  to  functions  defined  on  the  unit  sphere 
since  the  quadratic  form  grAg  is  a function  on  the  unit  sphere.  We  will  denote 
this  generalized  trace  operation  with  “gentr”.  For  functions  /(g),  with  antipodal 
symmetry  on  the  unit  sphere,  this  operation  is  given  by 


4.3  Theory 


(4.1) 


(4.2) 


(4.3) 


4.3.2  Mean  Diffusivity 
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The  mean  diffusivity  index  has  been  proposed  as  an  orientation  independent 
measure  of  diffusion,  in  the  context  of  traditional  (rank-2)  DTI  [57].  According 
to  its  original  definition,  it  is  given  by  the  mean  of  the  eigenvalues  of  the  rank-2 
diffusion  tensor,  which  is  just  the  trace  of  the  tensor  divided  by  3.  Since  diffusivity 
along  a direction  g,  assumed  by  traditional  DTI,  is  just  the  quadratic  form  of  the 
rank-2  diffusion  tensor,  D,  along  g (see  Eq.  2.32),  it  is  clear  from  Eq.  4.2  that  the 
original  definition  of  mean  diffusivity  is  just 

(D)  = Y [ D(g)dg  ■ (4.4) 

This  shows  that  the  mean  diffusivity  index  is  not  just  the  mean  value  of  the  three 
eigenvalues  or  that  of  the  diffusivities  along  three  orthogonal  directions,  but  it  is 
also  the  mean  value  of  diffusivities  along  all  directions  implied  by  rank-2  DTI.  Note 
that  comparison  of  Eqs.  4.3  and  4.4  implies  that  mean  value  of  a function  defined 
on  the  unit  sphere  is  just  one-third  of  its  generalized  trace: 

(/(g))  = ^gentr(/(g))  . (4.5) 

When  an  arbitrary  rank -l  Cartesian  tensor  is  used,  the  diffusivities  are  given 
by  Eq.  3.8.  As  a result  the  corresponding  mean  diffusivity  value  becomes 

1 3 3 3 

(0)  = 2;EE-Ed»* 

i i=l  42  — 1 ii— 1 

Here,  gik  is  the  4-th  component  of  the  vector  g as  given  in  Eq.  3.5.  The  integral 
in  the  above  expression  can  be  evaluated  analytically.  The  resulting  (D)  values  for 
ranks  up  to  6 are  provided  in  Table  4-1. 

Note  that  these  expressions  can  also  be  derived  by  incrementally  using  the 
expressions  relating  the  components  of  a higher  rank  tensor  to  those  of  lower  rank 
tensors  as  presented  in  Table  3-2.  This  is  done  by  using  the  expressions  in  this 
table  from  the  bottom  towards  the  top.  For  example,  the  diagonal  elements  of  the 
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Table  4 1:  Mean  diffusivity  values  in  terms  of  the  components  of  the  higher  rank 
tensors  through  rank  6. 


Rank 

(D) 

0 

D 

2 

| {Dxx  + Dyy  + Dzz) 

4 

^(-^xxxx  “1“  -^yyyy  Dzzzz  2 (-DXxyy  -t^xxzz  “I”  -Oyyzz)) 

6 

^(-C^XXXXXX  “1“  Dyyyyyy  + -C^ZZZZZZ 

""^(-Dxxxxyy  + -Dxxxxzz  H"  -Dyyyyxx  -Dyyyyzz  "t"  Dzzzzxx  + DZZZZy y) 

“hGZ^xxyyzz  ) 

rank-2  tensor  in  terms  of  the  components  of  the  rank-4  tensor  are  given  by  [77] : 

3 

^xx  - (9-f^xxxx  T 8DXXyy  T 8DXXZZ  Dyyyy  D ZZZ/  2 / \yZZ  ) 

Dyy  — ^(9-Dyyyy  T 8 Dxxyy  + 8 Dyyzz  — Dxxxx  ~ Dzzzz  — 2 Dxxzz)  (4.7) 

Dzz  = —(9  -Dzzzz  + 8DXXZZ  + 8Dyyzz  — Dxxxx  — 77yyyy  — 2Dxxyy) 


One-third  of  the  sum  of  these  three  expressions  is  just  the  rank-0  tensor  according 
to  the  first  line  of  this  table,  i.e., 


D — ~{DXX  + Dyy  + Dzz 


(4.8) 
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Note  that  this  is  the  expression  with  the  (D)  calculated  using  Eq.  4.6,  as  presented 
in  Table  4-1.  Therefore,  the  mean  diffusivity  value  is  just  the  diffusion  coefficient 
that  can  be  calculated  from  the  fitting  of  the  ffARDI  data  to  the  isotropic  Stejskal- 
Tanner  equation.  Furthermore,  since  the  same  result  can  be  obtained  by  starting 
from  any  of  the  rank-/  tensor,  mean  diffusivity  is  a model-independent  measure  of 
diffusivity. 

4.3.3  Anisotropy  as  the  Variance  of  the  Normalized  Diffusivities 

Despite  the  inflation  in  the  number  of  anisotropy  indices  already  proposed,  two 
of  them  (FA  and  RA)  have  been  the  ones  most  widely  used  [57].  These  measures 
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can  be  expressed  in  the  following  new  forms: 


fa = v K3  ” ■ and  ra = ^ trace(R2)  ■ 1 

Here,  R is  the  “normalized”  and  unitless  rank-2  diffusion  tensor, 


(4.9) 


R = 


D 


(4.10) 


trace(D) 

Therefore,  FA  and  RA  are  simply  functions  of  the  trace  of  the  square  of  this  tensor. 
In  terms  of  the  components  of  D,  this  quantity  is  given  by 


trace(R2)=^yj(BL  + + Dl  + 2 (D%  + D l + D%)).  (4.11) 

When  D is  diagonalized,  this  quantity  is  related  to  the  variance  of  the  eigenvalues 
of  the  normalized  rank-2  tensor.  Unlike  the  case  of  mean  diffusivity,  this  expression 
is  not  model  independent,  i.e.,  trace(R2)  is  not  equal  to  the  average  of  the  square 
of  the  diffusivities  along  all  directions.  In  this  work,  we  propose  to  use  the  variance 
of  the  normalized  diffusivities,  where  normalization  is  done  in  a similar  manner  to 
that  in  RA  and  FA,  and  show  how  this  measure  can  be  easily  generalized  to  higher 
rank  tensors  and  to  functions  defined  on  the  unit  sphere. 

We  start  by  defining  the  normalized  diffusivity  function  as 


Dn{  g) 


D(g) 

gentr(D(g))  ' 


(4.12) 


This  step  can  be  thought  of  as  the  generalization  of  the  transition  from  D to 
R in  Eq.  4.10.  Once  this  is  done,  instead  of  trace(R2),  we  can  use  the  quantity 
gentr(ZlAf(g)2).  When  a rank-/  tensor  model  is  used,  this  quantity  can  be  shown  to 
be  given  by 


gentr(Tbv(g)2) 
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where  Ni  = (l  + 1)(/  + 2)/2  is  the  number  of  unique  elements  of  the  rank-/  tensor, 
/rfcl  = l\/nx\ny\nz\  is  the  multiplicity  of  the  Aq-th  unique  element  of  the  diffusion 
tensor  Dkl,  and  gk1(pi)  is  the  component  of  the  unit  vector  specified  by  the  pi-th 
index  of  the  Aq-th  unique  element  of  the  diffusion  tensor  [77].  Note  that  in  the 
expression  for  /xkl,  nx,  ny  and  nz  are  respectively  the  number  of  x,  y and  2 indices 
in  the  full  sequence  of  subscripts  defining  the  component  of  the  tensor.  Definitions 
of  /ifc2  and  Dki  follow  the  same  lines  for  the  k2- th  component  of  the  tensor.  The 
integrals  in  Eq.  4.13  can  be  evaluated  analytically,  and  the  resulting  expressions 
are  listed  in  Table  4-2  for  tensors  up  to  rank-6. 

Having  evaluated  the  generalized  trace  of  the  square  of  the  normalized 
diffusivities  for  a rank -/  diffusion  tensor,  we  can  now  introduce  the  variance  of  the 
normalized  diffusivites  as  a measure  of  anisotropy.  This  is  done  by  using  Eq.  4.5: 


Because  this  is  the  variance  of  a normalized  distribution  of  diffusivities,  it  is 
unitless  and  has  no  diffusion  weighting. 

Variance  of  the  normalized  diffusivity  takes  its  minimum  value  of  0 only  when 
diffusivities  along  all  directions  are  equal.  When  a rank-/  tensor  model  is  used,  this 
constant  diffusivity  profile  is  achieved  when  all  terms  except  / = 0 in  its  irreducible 
representation  (Laplace  series)  are  0 [64].  Using  this,  we  have  listed  in  Table  4.3 
the  conditions  under  which  constant  diffusivity  profile  is  achieved  for  Cartesian 
tensors  up  to  rank-6.  Therefore,  regardless  of  whether  one  is  using  a rank-/  tensor 
model  or  a function  which  can  be  represented  in  general  with  a rank-00  tensor,  the 
theoretical  possible  value  for  minimum  variance  is  0 as  expected. 

Unlike  the  case  of  minimum  value,  the  supremum  value  of  the  variance  is  not 
the  same  in  different  rank  tensor  models.  It  is  possible  to  show  that  under  the 
condition  of  positive  semidefiniteness  of  the  tensor,  this  value  is  achieved  when 


V = variance^  (g))  = (DN( g)2)  - (DN{ g))2 


(4.14) 
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Table  4 2:  Generalized  trace  of  the  square  of  the  normalized  diffusivity  profiles  in 
terms  of  the  components  of  the  higher  rank  tensors  through  rank-6. 


Rank 

gentr  (D%) 

0 

~ I — 

3 

2 

45(D)2  ( 3(-DXx  + £yy  + ^zz)  + 4(^xy  + D'xz  + ^yz)  + 2(-Dxx.Dyy  + DXXDZZ  + DyyDzz)  ) 

4 

945(D)5  ( 35(HXxxx  + Dyyyy  + Dzzzz)  + 108(£>xxyy  + Aoczz  + Dyyzz) 

-(- 1 2 ( Dxxyy Dzzzz  4~  DXXzz  Dyyyy  “I"  Dyyzz  ^xxxx)  4~  60(DXXxxDXxyy 
~^DxxxxDxx  zz  4~  Dyyyy  DXXyy  “f*  Dyyyy  DyyZZ  4~  ^ZZZZ  ^XXZZ  4"  DzzzzDyyZZ) 

~\~72(DXXZZDyyZZ  “h  DXXyyDXXZZ  4"  D xxyy  Dy yZZ  ) 

+6(-DXxxx^yyyy  + DXXXXDZZZZ  4-  DyyyyDzzzz)  4-  144(Dxxyz  -f-  Dyyxz  4-  Dzzxy) 

+ 80(£>xXxy  -1-  Dxxxz  H-  Dyyyx  + Dyyyz  4"  D^  4"  D^y  ) 

+ 96(Z)XXXyi)yyyX  4“  D XXXy  D ZZXy  + D XXXZ  D yyXZ  4"  ^XXXZ^ZZZX  4"  D yyyXD  ZZXy 
+ DyyyZDZZZy  4"  DyyyzDXXyz  4~  D ZZZX  D yyXZ  4"  DZZZy  DXXyZ  ) ) 

6 

9009(D)2  ( ^31(HXXXXXX  + Dyyyyyy  + Dzzzzzz)  + 4860.Dxxyyzz 

4~10(DXXXXXX  Dyyyyyy  4“  DXXXXXX  DZZZZZZ  4"  D yyyyyy  D ZZZZZZ) 

+1575(Z?xxxxyy  + Dxxxxzz  + Dyyyyxx  + Dyyyyzz  + Dzzzzxx  + Dzzzzyy ) 

-f"630(Z)XXXXXXZ)XXXXyy  -|-  DXXXXXXDXXXXZZ  -h  Dyyyyyy  Dyyyyxx  4"  D yyyyyy  D yyyyzz 
~^DzzzzzzDzzzzxx  + DZZZiZZZDZZZZyy)  210(Z)XXXXXXZ)yyyyXX  4"  ^XXXXXX^ZZZZXX 

Dyyyyyy  DXXXXyy  + D yyyyyy  D ZZZZyy  + -^ZZZZZZ  ^XXXXZZ  4"  ^ZZZZZZ  D yyyyZZ  ) 
4"30(-DXXXXXXZ)yyyyzz  4"  ^XXXXXX  D ZZZZyy  4~  D yyyyyy  D xxxxzz  Dyyyyyy  D zzzzxx 

4~  DZZZZZZ  D XXXXyy  4-  DZZZZZZDyyyyXX)  + 2100  ( -^XXXXy  Z 4"  Dyyyyxz  4"  ^ZZZZX  y) 

4”420Z)XXyyZZ  (Z)xxxxxx  4”  Dyyyyyy  4“  -^ZZZZZz) 

4“1050(Z)XXXXyyZ)XXXXZZ  4“  DyyyyXXDyyyyZZ  D ZZZZXXD ZZZZ yy  ) 

4"450(Z)XXXXy  yZ)yyyyzZ  4"  ^XXXXZZ  D yyyyXX  D XXXXyy  D ZZZZXX  4“  D yyyyZZ  D ZZZZXX 

+ Dxxxxzz  -Dzzzzyy  4~  dyyyyxx  DZZZZyy  ) 

4"2700.DXXyyZZ  {DXXXXyy  4”  DxXXXZZ  4“  DyyyyXX  4”  DyyyyZZ  4“  D ZZZZXX  4”  D zzzz yy  ) 
4"270(DXXXXZZDyyyyZZ  4"  D XXXXyy  D ZZZZyy  4~  D yyyyXX  D ZZZZXX) 

4“ 2250  (DXXXXyy  Dyyyyxx  4"  D XXXXZZ  D zzzzxx  4”  DyyyyZZ  D ZZZZyy  ) 

+7o6(Z?xxxxxy  + Dxxxxxz  + Dyyyyyx  + Dyyyyyz  + Dzzzzzx  + Dzzzzzy) 

4-1680(DXXXXX2DXXXyyZ  4“  DXXXXXZDXXXZZZ  4~  DXXXyyyDyyyyyX  4“  DyyyXXZ  Dyyyyy  Z 

4"  Dyyyyyx  DyyyZZX  4"  Dyyyyy z DyyyZZZ  4"  D XXXXXy  D XXXyyy  D XXXXXy  D XXXZZy 

4”DZZZZZXDxXXZZZ  “1“  DZZZZZXDZZZ yyx  4“  DZZZZZy  DyyyZZZ  4“  DZZZZZy  DZZZXXy  ) 

4~360  (Dxxxxxz  Dyyyyxz  4“  Dyyyyyz  DXXXXyZ  4“  Dyyyyyx  D ZZZZXy  4"  D XXXXXy  Dyyyyy x 
4”  DxXXXXZ  DZZZZZX  4"  DZZZZZXDyyyyXZ  4"  D ZZZZZy  DXXXXyZ  4“  D yyyyyZ  D ZZZZZy 

4-DXxxxxyDZZZZXy)  4-  2000(Dxxxyyy  4-  Dxxxzzz  4-  Dyyyzzz) 

4-720(DXXXZZzDyyyyXZ  4"  DXXX  ZZy  Dyyyyy  X 4"  D XXXXyZ  D yyyZZZ  4"  Dy  yy  yy  Z D ZZZXXy 
4"  Dxxxxxz  DZZZyyX  4“  DXXxyyyDZZZZXy  4“  DXXxxxy  DyyyZZX  4“  DXXxyyzDzzzzzx 
4"  DyyyXXZ  DZZZZZy  ) 4“  43 20  ( DXXXZZy  Dyyyzzx  4“  D XXXyyZ  D ZZZyyX  4"  Dy  yy  XXZ  D ZZZXXy  ) 

+3600(.Dxxxyyz  + Dxxxzzy  + Dyyyxxz  + Dyyyzzx  + Dzzzxxy  + Dzzzyyx 

4"  DxXXXyZ  DyyyXXZ  4"  Dxxxyyz  Dyyyyxz  4"  DxXXXyZ  Dzzzxxy  4"  DyyyyXZDZZZyyX 
4”  Dxxxzzy  DZZZZXy  4-  DyyyZZXDZZZZXy)  4"  2400  (Dxxxyyy  Dxxxzzy  4"  DxxxyyZDxxxzzz 

4"DXXXyyy  DyyyZZX  4"  Dyyyxxz  DyyyZZZ  4"  D yyyZZZ  D ZZZXXy  4"  DXXXZZZ  DZZZy yX  ) ) 

the  tensor  is  given  by  a pure  outer  product  of  the  same  / vectors,  i.e.,  when  the 
components  of  the  tensor  are  given  by 


— D gilgi2  ■ ■ ■ git  . 


(4.15) 
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Table  4-3:  The  conditions  on  the  components  of  the  diffusion  tensors  that  yield 
isotropic  diffusivity  profiles. 


Rank 

Conditions 

0 

D is  finite. 

2 

Acx  — Dyy  = Dzz,  all  other  components  are  0. 

4 

CXxxx  = ^yyyy  = ^zzzz  — 3ZlXxyy  = 3 -DXxzz  = 3DyyZZ  , 

all  other  components  are  0. 

6 

DxXXXXX  -^yyyyyy  ^ZZZZZZ  5-DXXXXyy  — 5 -DXXXXZZ  — 5DyyyyXX 

5-DyyyyZZ  5-Dzzzzxx  5 -DZZZZyy  1 5 £^XXyyZZ  ? 

all  other  components  are  0. 

Here,  g'  is  the  unit  vector  specifying  the  direction  of  greatest  diffusion  coefficient 
where  D is  this  maximal  diffusivity.  Note  that  although  a real  generalized  diffusion 
tensor  may  come  arbitrarily  close  to  this,  it  can  never  reach  this  form,  since  it 
would  imply  0 diffusivities  along  directions  perpendicular  to  g'  as  can  be  seen  from 
Eq.  3.8.  Because  the  value  of  0 for  diffusivities  is  nonphysical,  we  refer  to  the 
variance  associated  with  the  tensor  given  in  Eq.  4.15  as  the  supremum  rather  than 
the  maximum  value.  It  is  possible  to  show  after  some  algebra  that  this  supremum 
value  corresponding  to  a rank -l  tensor  is  given  by 

l2 

sup  variance(L»N(g))  = + ^ . (4.16; 

This  result  indicates  that  the  supremum  value  is  dependent  on  the  rank  of  the 
tensor  model  selected.  This  is  a surprising  result  that  implies  that  there  is  an 
intrinsic  limit  to  the  anisotropy  that  can  be  quantified  with  a lower  rank  tensor 
model.  Moreover  for  functions  whose  domains  are  unit  hemisphere,  in  general,  the 
rank  of  the  tensor  model  required  for  an  exact  representation  of  these  functions 
is  infinity.  It  is  clear  from  Eq.  4.16  that  as  l — > oo,  the  supremum  value  goes  to 
infinity.  Therefore,  a reasonable  choice  of  function  that  will  quantify  anisotropy  in 
terms  of  the  variance  of  the  diffusivities  will  be  a monotonic  function  that  will  map 
the  interval  [0,  oo)  to  [0, 1).  Among  many  functions  that  obey  this  description,  we 
choose  to  use  an  expression  that  is  most  suitable  with  typical  datasets.  Upon  some 
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Figure  4-2:  The  continuous  curve  shows  the  GA  values  as  a function  of  the  vari- 
ance of  the  normalized  diffusivities  while  the  dashed  curve  shows  its  derivative. 
The  vertical  lines  indicate  when  the  supremum  values  of  the  variance  (hence  GA) 
are  achieved  for  tensors  of  ranks  2,  4 and  6 (from  left  to  right). 


scaling  we  arrive  at  the  definition  for  the  generalized  anisotropy: 

GA  = 1 - 


1 + (250V)dy)  ’ 


(4.17) 


where  the  exponent  e(V)  is  defined  as 


e(F)  = 1 + 


1 


1 + 50001/ 


(4.18) 


The  logarithmic  plot  of  GA  as  a function  of  the  variance  of  the  normalized  diffu- 
sivities is  shown  in  Fig.  4 2.  The  vertical  lines  in  this  plot  show  the  locations  of 
the  maximal  anisotropy  cases  that  can  be  quantified  by  tensors  of  ranks  2,  4 and  6 
respectively.  These  supremum  values  for  GA  are:  .957,  .980  and  .987. 

The  form  of  the  expression  for  GA  as  given  in  Eq.  4.17  is  not  arbitrary,  and 
takes  into  account  the  sensitivity  of  the  calculated  values  to  the  variations  in  the 
variance.  The  behavior  of  this  function  can  be  best  understood  by  examining  the 
derivative  of  this  expression  with  respect  to  V.  This  plot  is  also  presented  in  Fig. 

4 2 using  dashed  lines.  It  is  apparent  from  this  figure  that  the  formulation  of 
anisotropy  as  given  in  Eq.  4.17  yields  low  contrast  among  voxels  that  are  already 
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anisotropic  such  as  those  in  white-matter.  This  behavior  is  in  accordance  with  the 
previously  proposed  anisotropy  indices.  However,  unlike  the  previously  introduced 
indices,  GA  also  has  low  contrast  among  voxels  with  very  low  anisotropy  values 
such  as  those  in  free  water.  As  a result,  the  intensity  differences  in  GA  values  is 
concentrated  in  voxels  within  gray-matter  and  the  transition  from  gray-matter  to 
white-matter  while  retaining  high  intensity  in  white-matter.  If  one  is  interested  in 
changing  the  contrast  according  to  the  values  in  a different  kind  of  dataset,  this  can 
be  accomplished  with  ease  by  adjusting  the  constants  in  the  definition  of  GA. 

4.3.4  Isotropy  as  the  Entropy  of  the  Normalized  Diffusivity  Profile 

In  the  previous  section  we  showed  the  quantification  of  anisotropy  in  terms  of 
the  variance  of  the  diffusivity  values.  In  general  the  variance  in  a random  variable 
x is  given  by 


where  P(x)  is  the  probability  distribution  function  (PDF).  The  parametrization  of 
anisotropy  in  terms  of  the  variance  of  diffusivities  in  the  previous  section  effectively 
treated  the  diffusivity  as  a random  variable  and  assumed  a PDF  according  to  the 
frequencies  of  diffusivities  that  occur  in  the  diffusivity  profile.  An  alternative  way 
to  approach  the  same  problem  is  to  take  the  orientation  as  a random  variable  and 
use  Di\r( g)  as  the  PDF.  This  can  be  done  because  D^( g)  is  positive  definite  and 
integrates  to  1.  In  this  case,  the  quantity  that  we  refered  to  as  the  “variance” 
of  the  normalized  diffusivity  is  related  to  the  integral  of  the  square  of  a PDF. 

In  this  context,  it  is  not  clear  why  one  would  choose  to  do  this  instead  of  taking 
the  integral  of  any  other  power  of  the  PDF.  Perhaps  a more  appropriate  choice 
would  be  to  take  the  derivative  of  the  integral  with  respect  to  the  power  of  the 
distribution  and  set  this  power  equal  to  1.  We  will  refer  to  the  negative  of  this 


(4.19) 


quantity  as  a : 
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- - -1(Ud^^)L 

= ~Y  f DN(g)\nDN(g)dg  . (4.20) 

Therefore,  cr  is  a measure  of  how  the  integral  is  varying  when  the  power  of  the  PDF 
is  varied  around  1 (note  that  at  m = 1 the  integral  is  just  the  generalized  trace) 
and  is  equal  to  the  differential  entropy  of  the  distribution.  This  is  a well-known 
quantity  in  several  disciplines;  in  statistics,  it  quantifies  the  uniformness  of  the 
distribution  [78];  in  statistical  mechanics,  it  yields  the  uncertainty  level  in  a given 
phase  space  [79],  whereas  in  information  theory,  it  gives  the  information  content  of 
the  PDF  [78].  When  the  entropy  is  higher;  the  distribution  is  more  uniform,  the 
orientation  is  less  certain  and  we  have  less  orientational  information.  Therefore, 
when  anisotropy  is  viewed  as  a measure  of  certainty  or  information,  it  is  natural  to 
parametrize  it  in  terms  of  entropy. 

It  is  straightforward  to  show  that  for  an  isotropic  diffusivity  profile,  a will  take 
its  maximum  value  of  In  3.  On  the  other  hand  for  a rank -/  tensor  model,  as  the 
diffusivities  approach  the  form  given  in  Eq.  4.15,  entropy  approaches  its  infimum 
value  given  by 


inf  er  = 


/ 

1 + l 


Tin 


3 

1 + l 


(4.21) 


Obviously  for  an  arbitrary  function  which  can  in  general  be  represented  by  a rank- 
oo  tensor,  this  infimum  value  is  — oo.  Therefore,  a general  anisotropy  index  based 
on  entropy  has  to  be  a monotonically  decreasing  function  of  entropy  that  maps  the 
interval  (— oo,ln3]  to  [0,  1).  We  employ  a function  similar  to  that  in  Eq.  4.17,  and 
introduce  the  scaled  entropy  (SE)  index: 


SE  = 1- 


1 

1 + (60(ln3  - a))d>n3-<r)  ' 


(4.22) 


Here  the  exponent  function  e is  the  same  as  that  given  in  Eq.  4.18. 
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Figure  4 -3:  The  continuous  curve  shows  the  SE  values  as  a function  of  the  entropy 
of  the  normalized  diffusivities  while  the  dashed  curve  shows  its  derivative.  The  ver- 
tical lines  indicate  when  the  infimum  values  of  the  entropy  (hence  supremum  values 
of  SE)  are  achieved  for  tensors  of  ranks  2,  4 and  6 (from  right  to  left). 

The  exponential  plot  of  SE  and  the  derivative  of  SE  as  a function  of  the 
entropy  of  the  normalized  diffusivities  is  shown  in  Fig.  4 3.  The  vertical  lines 
in  this  plot  show  the  locations  of  the  maximal  anisotropy  profiles  that  can  be 
quantified  by  tensors  of  ranks  2,  4 and  6 respectively.  These  supremum  values  for 
SE  are:  .963,  .980  and  .985. 

4.3.5  Expressions  for  Arbitrary  Functions  Defined  on  the  Sphere 

Thus  far  we  have  treated  the  functions  defined  on  the  unit  sphere  as  the 
rank  — > oc  limit  of  the  higher  rank  tensors.  For  completeness,  we  present  some  of 
the  results  obtained  for  arbitrary  positive  definite  and  integrable  functions  below: 


Note  that  when  an  even  rank  tensor  is  used,  the  antipodal  symmetry  in  the  diffu- 
sivity  profile  is  automatically  achieved  [77] . Therefore,  in  our  previous  derivations 
it  was  sufficient  to  evaluate  the  integrals  on  the  hemisphere.  In  Eqs.  4.23-4.25,  we 


(4.23) 


(4.24) 


(4.25) 
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Figure  4 4:  The  calculated  scalars  shown  on  the  same  brain  slice.  Note  that  in 
the  case  of  rank-0,  there  is  no  orientational  variation  in  the  diffusivities,  so  the 
variance,  GA  and  SE  values  are  identically  0,  and  the  entropy  values  are  In  3 every- 
where in  the  a image. 


do  not  assume  the  antipodal  symmetry  of  the  function.  Therefore,  the  integrations 
are  on  the  whole  sphere,  S.  Once  the  (/(g)),  V and  a values  are  calculated,  the 
scaled  anisotropy  indices  GA  and  SE  can  be  calculated  using  Eqs.  4.17,  4.18  and 
4.22  as  before. 

4.4  Results  and  Discussion 

To  demonstrate  the  images  of  the  scalar  measures  derived,  we  have  acquired 
a series  of  diffusion  weighted  images  from  an  excised  rat  brain.  These  images  were 
fitted  to  the  generalized  Stejskal-Tanner  relation,  Eq.  3.14,  to  yield  the  components 
of  Cartesian  tensors  up  to  rank-6  using  methods  described  in  Chapter  3.  (D) 
values  were  evaluated  using  the  expressions  tabulated  in  Table  4-1.  Similarly,  exact 
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calculation  of  the  variance(£hv(g))  was  performed  by  using  the  expressions  in  Table 
4-2.  From  these,  GA  values  were  calculated  using  Eq.  4.17.  Calculation  of  the 
entropy  was  done  numerically  using  iterated  Gaussian  quadrature  employing  96 
transformation  points.  Then  SE  values  were  obtained  from  Eq.  4.22. 

Fig.  4 4 shows  the  images  calculated  from  this  rat  brain  at  a selected  slice.  As 
it  is  seen  from  the  first  column,  (D)  is  constant  regardless  of  the  rank  of  the  tensor 
model  being  used.  It  is  not  possible  to  use  rank-0  model  to  quantify  anisotropy. 
Although  images  from  different  rank  models  (rank  > 2)  look  similar  at  first  glance, 
it  is  possible  to  find  differences  between  images  from  different  rank  models.  GA 
and  SE  images  also  look  similar.  These  issues  will  be  discussed  below. 

4.4.1  Fiber  Crossings  and  Anisotropy 

Overcoming  the  inability  of  traditional  (rank-2)  DTI  to  resolve  multiple 
fiber  orientations  within  the  voxel  was  the  primary  motivation  to  generalize  the 
DTI  technique  to  employ  tensors  of  higher  ranks.  Therefore  an  appropriate  test 
for  the  application  of  generalized  DTI  should  be  the  illustration  of  changes  in 
the  anisotropy  values  when  models  with  different  rank  tensors  were  used.  More 
specifically,  one  can  expect  that  rank-2  DTI  will  yield  artificially  low  anisotropy 
values  in  pixels  with  complicated  fiber  structure.  To  test  this  hypothesis,  we  have 
performed  simulations  that  uses  the  exact  form  of  the  diffusional  signal  attenuation 
under  the  assumptions  that  the  fibers  are  perfect  cylinders,  and  water  molecules  are 
confined  to  the  interior  of  these  cylinders,  i.e. , the  membranes  are  not  permeable. 

In  this  case,  the  diffusion  weighted  MR  signal  attenuation  from  a cylinder  of  radius 
R and  length  L,  when  the  applied  diffusion  gradient  makes  an  angle  9 with  the 
direction  of  the  cylinder,  is  given  by  [80] 


S(  q) 
-So 


2KnmR2  (2'KqR)A  sin2(2 0)ot\m 
[( mrR/L )2  - (2tt qR cos 6)2}2 
[1  — (— l)T1cos(27rgLcos6')]  [J'rn(2'KqRs\.n9)f 

L 2 ialm  - fiirqR sin#)2]2  (a2km  - m2) 


x exp 


(4.26) 
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Table  4-4:  Simulated  GA  values  from  different  rank  tensor  models  for  1,  2 and  3 
fiber  orientations.  Also  included  are  the  mean  and  standard  deviations  of  the  GA 
values  that  were  calculated  from  any  region  of  the  simulated  image  presented  in 
Figure  4 1 containing  fibers. 


Rank- 2 

Rank- 4 

Rank-6 

1 fiber 

0.89037 

0.89036 

0.89035 

2 fibers 

0.56322 

0.63429 

0.63419 

3 fibers 

1.19  x 10"7 

0.19548 

0.19914 

Image  in  Fig.  4 1 

0.844  ± 0.101 

0.852  ± 0.081 

0.852  ± 0.082 

In  this  expression,  q = (27r)_17JG,  Jm  is  the  m-th  order  Bessel  function,  a^m  is 
the  k- th  solution  of  the  equation  J'm(ot ) = 0 with  the  convention  a10  = 0,  and 
A = tinodmo  + 2 [(1  — Sn0)  + (1  — <5mo)],  where  is  the  Kronecker  delta.  In  the 
presence  of  more  than  one  cylinder,  the  signal  attenuations  from  these  cylinders 
become  additive  if  the  signals  are  equal  when  no  diffusion  gradient  is  applied. 

We  have  evaluated  Eq.  4.26,  with  the  parameters:  L = 5 mm,  R = 5/rm, 

D = 2.0  x 10 _3mm2/s,  A = 17.8ms,  5 = 2.2ms,  b = 1500 s/mm2.  These  resemble 
our  MRI  experiment  on  the  excised  rat  brain.  Similar  to  that  in  Hagen  and 
Henkelman  [69],  we  have  terminated  the  infinite  series  at  n = 1000  and  fc,m  = 10. 
Just  as  in  our  HARDI  experiment,  the  gradient  directions  were  chosen  to  point 
toward  the  81  vertices  of  the  tessellations  of  an  icosahedron  on  a unit  hemisphere 
while  keeping  the  fiber  direction  fixed.  This  way,  a number  of  signal  attenuations 
are  computed.  Then  we  use  Eq.  3.14  to  compute  the  components  of  the  tensors  up 
to  rank-6.  From  these,  we  calculate  the  GA  values  for  1,  2 and  3 fiber  cases  when 
the  angle  between  any  two  distinct  fiber  direction  is  90°.  The  results  are  tabulated 
in  Table  4-4.  Using  the  same  idea,  we  have  constructed  a simulation  image  matrix 
of  256  x 256,  which  was  already  presented  in  Figure  4 1.  This  simulation  image 
contains  two  fiber  bundles  one  of  which  is  linear,  the  other  is  circular.  This  way,  a 
distribution  of  angles  is  obtained  in  different  pixels  of  the  image.  In  Table  4-4,  we 
also  give  the  mean  GA  values  and  their  standard  deviations  calculated  from  the 
section  of  this  image  containing  fibers. 
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It  is  clear  from  the  first  row  of  Table  4-4  that  in  a voxel  containing  a single 
fiber  orientation,  the  measured  anisotropy  is  almost  independent  of  the  rank  of 
the  model  being  used.  However,  the  situation  is  quite  different  when  there  are 
more  fiber  orientations.  When  there  are  two  fiber  orientations,  the  difference  in 
the  calculated  GA  values  is  significantly  higher  in  rank-4  and  rank-6  models  when 
compared  to  a rank-2  model.  This  change  is  even  more  dramatic  in  the  case  of 
three  fiber  orientations.  In  this  case,  the  GA  value  implied  by  the  rank-2  tensor 
is  so  small  that  one  can  easily  label  the  geometry  as  non-fibrous,  where  higher 
rank  tensors  make  the  anisotropy  apparent.  The  values  from  the  simulated  image 
shown  in  Figure  4-  1,  presented  in  the  last  line  of  Table  4-4,  are  also  consistent  with 
these  findings.  Just  like  one  would  specify  a region  of  interest  in  a real  dataset  to 
calculate  the  anisotropy  values  from  a white-matter  structure,  the  section  of  this 
simulation  image  that  contains  one  or  two  fibers  are  included  in  the  calculation  of 
the  GA  values.  It  should  be  noted  that  in  addition  to  the  elevation  of  the  mean 
values,  there  is  also  a decrease  in  the  standard  deviations  of  the  GA  values  when 
higher  rank  tensor  models  are  used.  This  is  because  as  the  rank  of  the  tensor 
model  is  increased,  the  anisotropy  values  in  the  region  with  crossing  fibers  approach 
GA  values  from  those  voxels  with  single  orientations,  decreasing  the  overall  spread 
in  the  values.  For  example,  as  can  be  seen  from  the  second  line  of  Table  4-4,  in 
the  case  of  two  fibers,  the  GA  value  rises  from  0.56  to  0.63  when  rank  of  the  tensor 
is  increased,  driving  the  GA  values  towards  the  GA  value  from  voxels  with  single 
fiber,  which  is  0.89.  When  we  look  at  just  the  section  with  crossing  fibers  in  Figure 
4-1,  we  see  that  the  enhancement  in  the  mean  of  the  anisotropy  values  is  larger 
than  7%,  which  is  quite  significant.  These  findings  can  also  be  seen  in  Figure  4 5, 
where  we  present  the  GA  images  from  rank-2  and  rank-6  tensor  fits.  Also  included 
in  this  figure  are  the  differences  in  variance  and  GA  values  from  these  two  tensor 


models. 
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GA(rank-2)  GA(rank-6)  AV  AGA 


Figure  4 5:  From  left  to  right:  GA  images  from  rank-2  and  rank-6  tensors,  the  in- 
crease in  the  variance  and  GA  values  when  the  rank  of  the  tensor  model  was  raised 
from  2 to  6.  All  images  are  calculated  for  the  simulated  system  shown  in  Fig.  4-1. 

Sensitivity  of  the  anisotropy  values  to  changing  ranks  provides  a potential  tool 
in  the  identification  of  voxels  with  heterogeneity  in  fiber  orientations.  However, 
one  should  keep  in  mind  that  anisotropy  indices  such  as  GA  and  SE  are  functions 
of  some  other  quantities  (variance  and  entropy)  and  the  changes  in  the  values  of 
these  indices  are  dependent  not  only  on  the  structural  complexity  within  the  voxel, 
but  also  on  how  these  anisotropy  measures  are  defined.  For  example,  since  the 
derivative  of  GA  has  a peak  around  the  variance  value  of  about  10-3,  small  changes 
in  the  voxels,  with  variance  values  in  this  range,  may  be  exaggerated  in  the  GA 
difference  maps  calculated  using  different  tensor  ranks.  For  this  reason,  one  may 
think  that  difference  of  the  variance  or  entropy  maps  (or  some  other  functions  of 
these)  may  be  a more  suitable  choice  when  one  is  interested  in  identifying  regions 
with  complex  fiber  architecture.  However,  note  that  when  variance  or  entropy  is 
used  for  this  purpose,  the  supremum  or  infimum  values  differ  with  changing  ranks 
because  as  mentioned  previously,  there  is  a limit  to  the  anisotropy  that  can  be 
quantified  by  a given  model.  Therefore,  one  may  get  highly  different  values  of 
variance  or  entropy  from  different  tensor  models  in  voxels  with  very  anisotropic 
structure.  As  a result,  in  this  case  one  should  be  interested  in  voxels  with  high 
variance  (entropy)  differences  but  not  so  high  anisotropy  values  because  the  limits 
of  both  variance  and  entropy  values  at  high  anisotropy  depend  on  the  rank. 
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Figure  4-6:  GA  values  from  rank-2  (left  column)  and  rank-6  (second  column) 
tensors  from  different  coronal  slices.  The  right  two  columns  show  the  difference 
between  the  variance  and  GA  values  when  these  two  tensor  models  were  used. 


Based  on  these  findings  from  simulations,  we  have  investigated  how  the  GA 
values  in  tissue  change  with  changing  tensor  ranks.  For  this  purpose,  we  show 
images  of  GA  values  from  rank-2  and  rank-6  tensors  in  Fig.  4-6  along  with  the 
difference  maps  of  the  variance  and  GA  values  when  these  two  tensor  models 
are  employed.  The  difference  images  are  scaled  so  that  the  maximum  intensity 
corresponds  to  a GA  difference  of  0.1  and  a variance  difference  of  0.001.  As  it 
is  clear  from  the  difference  maps  in  the  right  two  columns  of  Fig.  4-6,  there  are 
certain  regions  where  these  differences  are  significant.  In  the  first  two  rows,  these 
regions  are  within  the  pons,  and  in  the  bottom  row,  the  GA  values  in  the  brain 
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stem  are  enhanced.  These  regions  are  known  to  contain  structures  with  crossing 
fibers.  The  difference  images  present  significant  level  of  symmetry  implying  that 
the  anisotropy  and  variance  changes  are  most  likely  not  due  to  noise.  Rather, 
they  stem  from  local  structure  (such  as  crossing  fibers)  and  partial  volume  effects. 
Another  point  to  note  is  that,  as  it  was  noticed  in  simulations,  the  GA  images 
from  higher  rank  tensor  model  is  smoother  than  those  from  lower  rank  tensors. 
Although  the  diffusivity  profiles  from  individual  pixels  are  getting  sharper,  the 
resulting  anisotropy  maps  are  getting  smoother  as  a result  of  the  more  accurate 
quantification  of  anisotropy  when  the  rank  of  the  tensor  model  is  increased. 

4.4.2  Simple  Implementation  of  Entropy  Based  Anisotropy  for  Rank-2  Tensors 

The  underlying  idea  in  the  construction  of  FA  and  RA  indices  was  that 
anisotropy  is  a quantity  that  depends  upon  the  variance  of  the  normalized  eigenval- 
ues, which  are  the  normalized  diffusivities  along  the  principal  axes  of  the  diffusion 
ellipsoid.  We  showed  in  Eq.  4.9  that  these  indices  can  be  expressed  in  terms  of 
trace(R2).  Our  generalization  of  these  indices  included  the  extension  of  the  vari- 
ance to  diffusivities  along  all  directions,  and  generalization  of  the  trace  operator.  If 
we  look  at  those  rank-2  measures  from  the  perspective  of  our  generalization,  scalar 
measures  like  FA  and  RA  can  be  thought  of  as  a simplification  to  the  GA  values  for 
the  case  of  rank-2. 

Now,  we  consider  the  formulation  of  SE  index  and  ask  the  question:  Is  there 
a simplified  formulation  of  entropy  for  rank-2  tensors?  This  would  be  particularly 
useful  in  the  case  of  SE  index,  because  the  numerical  implementation  of  the 
integral  given  in  Eq.  4.20  is  a slow  computational  operation.  Since  the  diffusivity 
along  a direction  g assumed  by  the  rank-2  tensor  model  is  just  the  quadratic  form 
of  the  diffusion  tensor  as  given  in  Eq.  2.47,  the  normalized  diffusivity  is  just 


Dn(s)  = grRg 


(4.27) 


77 

where  R is  defined  in  Eq.  4.10.  Therefore,  entropy  of  a rank-2  diffusion  tensor 
becomes  (from  Eq.  4.20) 

a = gTRg  ln(gTRg)  dg  . (4.28) 

Now,  we  use  an  approximation  that  is  well-known  in  the  field  of  mathematical 
physics  that  establishes  the  relation  between  classical  and  quantum-mechanical 
entropy  [81]: 

y gTRln(R)gfl!g  • (4.29) 

Here,  the  In  operation  of  the  matrix  R is  defined  in  terms  of  the  Taylor  expansion 

oo 

lnR=  -^-(I-R)n  , (4.30) 

71=1  ^ 

where  I is  the  identity  matrix.  From  Eq.  4.1,  it  is  easy  to  see  that 

a ~ — trace(RlnR)  . (4.31) 

The  right  hand  side  of  the  above  equation  is  the  same  expression  as  the  “von 
Neumann  entropy”  of  a density  matrix  (that  will  be  denoted  by  <tvn),  well  known 
in  quantum  mechanics.  The  evaluation  of  the  infinite  summation  in  Taylor 
expansion  in  Eq.  4.30  can  be  avoided,  since  upon  diagonalization  of  R,  ctvn  is 
given  by 

3 

°vN  = A In  # , (4.32) 

2=1 

where  p,  are  the  eigenvalues  of  R,  or  eigenvalues  of  the  original  diffusion  tensor  D 
divided  by  the  trace  of  the  tensor. 

The  “approximation”  used  in  Eq.  4.29  is  an  unusual  one,  which  is  very 
accurate  for  reasonably  isotropic  tensors.  However,  as  the  tensor  approaches  the 
most  anisotropic  limit,  <tvn  deviates  from  a,  as  it  is  obvious  from  their  infimum 
values  of  0 for  ctvn  and  2/3  for  cr  in  the  case  of  rank-2  [81].  Note  that  a similar 
situation  occurs  in  the  comparison  of  trace  (R2)  with  a supremum  value  of  1 and 
gentr(H/v(g)2)  with  a supremum  value  of  3/5  for  rank-2.  Therefore  just  as  FA 
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Figure  4-7:  Dependence  of  the  FA,  GA  and  VA  (from  left  to  right)  values  on  the 
two  larger  eigenvalues  of  the  matrix  R (top  row),  and  images  of  these  scalar  mea- 
sures from  the  same  slice  of  the  rat  brain  (bottom  row). 


values  require  a scaling  different  from  that  of  GA  values,  scaling  of  ervN  to  an 
anisotropy  measure  will  be  different  from  that  of  a.  Following  a scaling  similar  to 
the  one  used  in  the  formulation  of  FA,  we  arrive  at  the  definition  for  the  anisotropy 
measure  that  we  have  called  visual  anisotropy  [82]: 

VA=f^T( 2S3TT-  (4-33) 

These  formulations  enable  one  to  visualize  the  behavior  of  these  indices  for 
given  eigenvalues  of  R.  If  we  order  the  eigenvalues  such  that  pi  > p2  > Pz , then 
this  condition,  along  with  trace(R)  = p\  + p2  + pj,  — 1,  and  the  positive  definiteness 
of  the  eigenvalues,  limits  the  allowed  values  of  p\  and  p2  to  the  triangular  zone  in 
the  pip2  plane  given  by  pi  - p2  > 0,  pi  + 2p2  > 1,  p\  + p2  < 1.  (For  a similar 
visualization  see  Alexander  et  al.  [83].)  These  figures  along  with  images  of  FA, 

GA,  VA  indices  are  shown  in  Fig.  4 7.  Note  that  the  difference  in  the  brightness  of 
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these  images  are  due  to  the  different  scalings  that  were  used  in  the  computations  of 
the  indices. 

4.4.3  Variance  or  Entropy  Based  Anisotropy? 

In  this  work,  we  have  presented  two  different  ways  to  quantify  anisotropy: 
defining  it  in  terms  of  variance  and  entropy.  We  showed  how  both  of  these  ap- 
proaches can  be  generalized  to  higher  rank  tensor  models  as  well  as  to  model- 
independent  functions  (rank-oo  tensors)  defined  on  the  sphere.  Although  these  two 
formulations  provide  similar  looking  images,  it  is  possible  to  show  that  they  differ 
when  one  tries  to  order  pixels  according  to  their  anisotropy.  In  other  words  when 
one  compares  the  anisotropy  values  from  two  pixels,  the  GA  values  may  suggest 
that  the  first  pixel  is  more  anisotropic  than  the  second  while  SE  index  may  suggest 
otherwise.  This  disagreement  stems  from  the  difference  in  the  meanings  attached 
to  anisotropy.  Variance  is  a measure  of  dispersion,  while  entropy  is  a measure  of 
(lack  of)  information  or  uncertainty.  The  choice  between  the  two  depends  on  the 
utilization  of  the  anisotropy  index.  For  example,  using  the  fall  of  anisotropy  as 
a termination  criterion  in  fiber-tract  mapping  algorithms  assumes  that  the  fiber 
direction  is  too  uncertain  in  these  regions,  implying  that  anisotropy  index  being 
used  is  a measure  of  uncertainty. 

Because  of  the  logarithm  in  the  definition  of  entropy,  it  is  advantageous  to  use 
GA  from  a practical  point  of  view.  As  we  showed  in  Eq.  4.17  and  Table  4-2,  it  is 
possible  to  analytically  evaluate  the  GA  values  for  higher-rank  tensors.  The  ana- 
lytical evaluation  of  integrals  in  the  calculation  of  entropy  is  not  straightforward. 
Therefore,  we  adopted  a numerical  approach.  As  a result,  calculation  of  GA  values 
is  exact  and  faster  when  compared  to  the  calculation  of  SE  values.  When  one  deals 
with  model  independent  functions,  it  will  probably  be  necessary  to  evaluate  the 
variance  numerically  as  well. 


4.4.4  Possible  Extensions 
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Defining  anisotropy  as  a measure  of  information  in  a PDF  makes  it  possible  to 
come  up  with  different  measures  which  are  known  from  information  theory.  As  an 
example  the  relative  entropy  is  a measure  of  closeness  of  two  distributions  [79].  If 
we  associate  two  distributions  with  two  normalized  positive  definite  functions  (/1Ar 
and  f2N)  on  the  unit  sphere,  then  relative  entropy  (also  known  as  I<L  divergence)  is 
given  by 

ff(/uv(g)ll/j»(g))  = Jj:  jf /.*(g)ln|^ rfg  ■ (4.34) 

In  the  literature,  there  are  also  distance  measures  that  define  how  far  the  two 
distributions  are  apart.  First  of  these  is  called  Kolmogorov  distance  given  by 

^(/uv(g)./2Jv(g))  = g-  f |/iw(g)  - /2jv(g)|  dg  , (4.35) 

J s 

and  fidelity  defined  as 

Wuv(g),/2jv(g))  = J-  j v7iJv(g)/2Jv(g)dg  • (4.36) 

J s 

A detailed  analysis  of  these  measures  as  well  as  their  quantum  analogs,  in  which 
case  a density  matrix  analogy  of  normalized  rank-2  diffusion  tensor  can  be  used, 
can  be  found  in  Nielsen  and  Chuang  [84], 

Instead  of  using  the  data  as  the  PDF,  one  can  also  construct  a PDF  from 
the  histogram  formed  from  the  data  values  which  was  intrinsically  assumed  in  the 
variance  based  measures.  Given  two  such  functions  (such  as  in  two  voxels),  /i(g) 
and  f2( g),  one  can  construct  the  individual  PDFs  p(fl)  and  p(f2),  as  well  as  the 
joint  PDF  from  the  two  dimensional  histogram  of  the  two  functions  p(fi,f2).  In 
this  case,  the  KL  divergence  of  the  joint  probability,  and  the  product  of  the  two 
PDFs  is  called  the  mutual  information: 


Mi(/1,/2)  = i/(p(/i,/2)ib(/1)P(/2)) 


(4.37) 
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where  H is  defined  in  Eq.  4.34.  These  scalar  measures  given  in  Eqs.  4.34-4.37  may 
be  found  useful  in  applications  that  require  a similarity  index  between  different 
tensors  or  functions  such  as  image  smoothing  and  registration. 

The  results  we  have  presented  for  the  scalar  measures  derived  in  this  paper 
used  the  diffusivity  profile.  However,  the  formulations  also  remain  valid  for  other 
applications  such  as  g-space  imaging  in  which  one  calculates  the  average  propagator 
for  water  molecules  if  an  orientation  distribution  function  (ODF)  is  constructed,  as 
in  Tuch  et  al.  [85],  by  projecting  the  average  propagator  radially  on  to  the  surface 
of  a sphere.  This  is  because  both  the  ODF  and  normalized  diffusivity  profiles  are 
defined  on  the  unit  sphere,  they  are  positive  definite  and  integrate  to  1.  Note  that 
our  analysis  started  with  a non-normalized  diffusivity  profile  and  normalized  it 
using  the  generalization  of  the  trace  operation.  Therefore,  the  formulations  are 
applicable  to  all  kinds  of  positive  valued  data  acquired  on  or  projected  on  to  the 
surface  of  a sphere. 

4.5  Summary  and  Conclusions 

We  have  presented  rotationally  invariant  indices  of  mean  diffusivity  ((D))  and 
anisotropy  based  on  variance  and  entropy  (GA  and  SE).  Our  formulations  can  be 
used  to  quantify  mean  diffusivity  and  anisotropy  from  higher  rank  diffusion  tensors 
as  well  as  from  arbitrary  functions  defined  on  the  unit  sphere. 

Our  results  indicate  that  in  the  quantification  of  ( D ) it  is  appropriate  to  use 
any  rank  diffusion  tensor.  However,  lower  rank  tensors  may  undesirably  suppress 
the  anisotropy  information  that  is  available  from  HARDI  experiments.  This  fact 
was  predicted  through  simulations  and  confirmed  with  experiments. 


CHAPTER  5 

FIBER  ORIENTATION  MAPPING  IN  COMPLICATED  GEOMETRIES 

5.1  Preface 

Generalized  diffusion  tensor  imaging  uses  tensors  of  ranks  higher  than  two  to 
model  the  angular  variations  in  the  diffusivities  measured  by  magnetic  resonance 
imaging  (MRI)  methods.  However,  a diffusivity  profile  alone  is  not  readily  capable 
of  producing  distinct  fiber  orientations.  In  this  chapter,  we  show  how  it  is  possible 
to  approximate  the  displacement  probability  profile  for  water  molecules  from  the 
higher  rank  diffusion  tensors  and  validate  the  technique  via  simulations  of  one,  two 
and  three  fiber  systems.  Finally,  we  present  fiber  orientation  results  from  images 
of  an  excised  rat  brain  and  spinal  cord.  We  also  present  how  the  orientation  maps 
can  be  estimated  directly  from  a Laplace  series  expansion  of  the  displacement 
probabilities.  Using  this  technique,  it  is  possible  to  generate  similar  orientation 
maps  in  a much  shorter  computational  time  frame. 

5.2  Introduction 

Water  molecules  in  fibrous  tissues  like  white-matter  diffuses  mainly  along 
directions  parallel  to  the  fibers.  This  fact  can  be  used  to  infer  information  about 
connectivity  of  different  regions  in  the  brain.  The  most  common  approach  taken 
to  this  end,  DT-MRI,  employs  a second  order  diffusion  tensor  [51]  in  the  equation 
describing  the  evolution  of  magnetization  (Bloch- Torrey  equation)  [17].  We  have 
shown  in  Chapter  2,  how  one  can  estimate  the  fiber  orientation  starting  from  the 
diffusion  attenuated  multidirectional  signal.  Although  this  approach  approximates 
the  displacement  profile  to  an  oriented  Gaussian  function,  which  may  be  an 
oversimplification  in  most  cases,  the  generated  orientation  maps  using  DT-MRI 
have  been  found  quite  accurate  when  there  is  a single  orientation  in  the  region  of 
interest.  We  have  seen  in  Eqs.  2.50-2.52  that  DT-MRI  achieves  this  by  assuming 
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monoexponential  attenuation  along  each  orientation,  where  the  decay  rate  along  a 
particular  direction  is  determined  by  a diffusion  coefficient  which  is  the  quadratic 
form  of  the  tensor  along  that  direction.  However,  the  oriented  Gaussian  functions 
implied  by  DT-MRI  can  not  accommodate  more  than  one  fiber  directions.  This 
reduces  the  reliability  of  orientation  maps  since  one  will  likely  get  orientational 
heterogeneity  in  many  voxels  due  to  several  factors  such  as  fiber  crossings  and 
partial  volume  effects. 

Our  generalization  of  DT-MRI  is  capable  of  producing  diffusivity  profiles 
that  are  more  complicated  than  the  ones  that  can  be  constructed  from  rank-2 
tensors.  Although  this  scheme  yielded  the  diffusivities  along  all  directions,  the  fiber 
directions  in  general  are  not  directly  achievable  from  a higher  rank  diffusion  tensor. 
In  Section  5.3  we  attempt  to  generate  orientation  maps  implied  by  the  generalized 
(rank  > 2)  tensors  using  a numerical  approach.  In  Section  5.4,  we  employ  a more 
direct  approach  and  attempt  to  calculate  the  orientation  maps  from  the  signal 
attenuation  values  without  going  through  the  Cartesian  tensor  fitting  process. 


Similar  to  the  case  of  traditional  (rank  2)  DTI,  in  this  section  we  assume 
that  the  signal  decay  is  monoexponential  along  any  direction  but  with  a diffusion 
rate  specified  by  Eq.  3.8  instead  of  Eq.  2.47.  In  this  case,  the  average  propagator 
becomes  (using  Eq.  2.51) 


for  l > 2.  Therefore,  we  have  adopted  a numerical  scheme.  In  this  approach  the 
signal  values  are  synthetically  generated  on  a three  dimensional  Cartesian  lattice  in 
g-space.  This  is  done  by  using  the  monoexponentiality  assumption 


5.3  Orientations  From  Generalized  DTI 


(5.1) 


Analytical  evaluation  of  the  integral  in  the  above  expression  is  a formidable  task 


= exp  (-47tV  (A  - 6/3)  D( g))  , 


(5.2) 
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Figure  5-1:  Simulations  of  three  voxels  containing  one,  two  and  three  fiber  ori- 
entations (from  top  to  bottom).  The  left  column  shows  the  fiber  system  being 
simulated.  The  second  column  shows  the  ellipsoids  implied  by  traditional  DTI. 
ADC  profiles  are  depicted  in  the  middle  column.  Final  two  columns  demonstrate 
the  results  from  rank-8  generalized  DTI. 

where  D{ g)  is  given  by  Eq.  3.8.  Then  a fast  Fourier  transform  (FFT)  algorithm 
[86]  was  employed  to  generate  the  average  propagator.  The  isosurfaces  of  these 
propagators  are  calculated  and  visualized. 

5.3.1  Simulation  Results 

The  above  procedure  was  applied  to  the  simulations  of  the  MR  signal  in 
cylindirical  geometry  as  presented  in  Eq.  4.26.  One,  two  and  three  fiber  systems 
were  considered.  The  synthetically  generated  g-space  values  were  calculated  on  a 
64  x 64  x 64  grid  such  that  the  largest  g-value  corresponds  to  a 6-value  of  60000 
s/mm2.  The  results  of  this  scheme  from  these  systems  is  presented  in  Figure  5 1 
along  with  the  ADC  profiles  and  diffusion  ellipsoids  generated  using  rank-2  DTI. 

It  is  clear  that  rank-2  DTI  performs  well  in  one  fiber  system,  but  it  is  not  able 
to  resolve  more  than  one  fiber  directions  in  complicated  systems.  Similarly,  ADC 
profiles  also  give  the  correct  orientational  information  when  there  is  only  one  fiber 
direction.  Although  the  ADC  profiles  get  more  complicated  as  the  number  of  fiber 
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Figure  5 2:  Calculation  of  fiber  probability  surfaces  from  the  water  displacement 
equiprobability  surfaces.  The  second  row  shows  that  an  ellipse  will  turn  into  a 
‘peanut’  shaped  object  under  this  transformation. 


orientations  is  increased,  the  peaks  of  the  ADC  profiles  no  longer  yield  the  distinct 
fiber  orientations  in  systems  with  orientational  heterogeneity.  The  last  two  columns 
in  this  image  depict  the  results  obtained  using  generalized  DTI  employing  rank-8 
tensors.  The  first  of  these  columns  show  the  isosurfaces  of  the  average  propagator 
calculated  from 

P(AR)  = PQ  . (5.3) 

Clearly,  the  peaks  of  this  surface  give  accurate  information  on  the  number  and 
orientations  of  the  distinct  fibers.  The  last  column  show  the  same  surface  after  a 
sharpening  transformation  that  we  have  implemented  to  make  the  results  easier  to 
interpret.  This  transformation  involves  the  “subtraction”  of  the  largest  sphere  that 
fits  into  the  isosurface.  This  can  be  thought  of  as  the  analog  of  keeping  only  the 
principal  eigenvector  in  rank-2  DTI.  This  transformation  is  depicted  in  Figure  5-2 
for  one  dimensional  surfaces. 

Note  that  in  the  above  description,  one  has  the  option  of  setting  the  proba- 
bility to  a certain  value  (P0)  when  constructing  the  isosurfaces.  We  have  found 
that  the  locations  of  the  peaks  of  the  isosurfaces  are  not  affected  by  the  choice  of 
P0  over  a large  range.  Figure  5 3 shows  how  the  surfaces  evolve  with  increasing 
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Figure  5 3:  Sharpened  probability  isosurfaces.  These  equiprobability  surfaces  were 
constructed  by  setting  the  probability  value  (P0)  from  1 x 10~4  to  7 x 10-4  in 
equal  steps  (from  left  to  right).  Top  row  shows  these  surfaces  when  there  is  a single 
orientation,  where  the  bottom  row  shows  the  same  surfaces  from  a voxel  with  two 
distinct  orientations. 


probability  values.  Note  that  increasing  the  probability  value  amounts  to  con- 
structing the  surface  closer  to  the  origin,  similar  to  the  short  time  behavior  of 
displacements.  Therefore,  the  constructed  surfaces  are  getting  smoother  as  the  P0 
value  is  increased.  Similarly,  one  may  expect  the  sharpness  of  these  surfaces  to  be 
dependent  upon  the  acquisition  parameters,  particularly  diffusion  time  and  6-value. 

Next  we  apply  the  same  techniques  to  the  synthetic  image  that  was  presented 
in  Figure  4-1.  The  image  matrix  was  reduced  to  16  x 16  for  clarity.  A rank-6  tensor 
model  was  used.  The  resulting  orientation  surfaces  are  shown  in  Figure  5 4.  In  the 
construction  of  this  image  P0  was  set  to  1 x 10-4.  Note  that  as  the  angle  between 
the  distinct  fibers  decrease,  the  two  peaks  tend  to  merge.  This  indicates  that  it 
may  not  possible  to  resolve  small  angular  deviations. 

5.3.2  Imaging  Results 

We  have  applied  the  procedure  detailed  above  to  multidirectional  diffusion 
weighted  imaging  data  collected  from  excised  rat  brain.  Similar  to  the  analysis 
performed  in  the  simulations,  signal  values  on  a Cartesian  grid  whose  points 
correspond  to  different  gradient  values  (q-space)  were  calculated  assuming  that 
signal  is  exponentially  decaying  along  each  radial  line.  This  way,  a 32  x 32  x 32 
g-space  image  covering  a region  up  to  6 = 80000s/mm2  was  constructed,  whose 
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Figure  5-4:  Isosurfaces  of  the  water  displacement  probability  functions  after  the 
sharpening  transformation.  The  geometry  is  the  subsampled  version  of  the  system 
presented  in  Figure  4-1. 

Fourier  transform  yielded  the  displacement  probabilities.  The  isosurfaces  of  these 
displacement  probabilities  were  transformed  using  the  transformation  depicted  in 
Fig.  5-2.  These  surfaces  were  examined  at  selected  slices  and  in  several  regions  of 
interest  (ROIs). 

Fig.  5 5 shows  the  fiber  structure  in  two  different  ROIs  in  a rat  brain.  The 
image  in  the  middle  is  the  fractional  anisotropy  [57]  map  constructed  from  a rank-2 
diffusion  tensor,  and  shows  the  location  of  the  selected  ROIs  in  the  brain.  The 
image  on  the  left  is  the  transformed  equiprobability  surfaces  in  the  ROI  on  the 
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Figure  5-5:  Isosurfaces  of  the  water  displacement  probability  functions  after  the 
sharpening  transformation.  The  middle  image  shows  two  ROIs  drawn  on  an  FA 
map.  The  images  on  the  sides  depict  the  orientation  isosurfaces  calculated  from 
these  (3.15mm)2  regions. 


left  of  the  FA  map,  and  the  third  image  shows  the  fiber  structure  for  the  ROI  on 
the  right  of  the  FA  map.  It  is  apparent  from  the  probability  maps  that  using  the 
generalized  diffusion  tensor  imaging,  one  is  able  to  produce  probability  surfaces 
which  are  clearly  non-Gaussian.  Note  that  in  the  major  white-matter  pathways, 
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Figure  5-6:  Isosurfaces  of  the  water  displacement  probability  functions  after  the 
sharpening  transformation.  The  image  on  the  upper  right  corner  is  the  image  with 
no  diffusion  weighting.  The  orientation  surfaces  calculated  from  the  selected  ROI 
(a  1.44mm  by  2.22 mm  region)  on  this  image  is  shown  on  the  left. 


the  generated  surfaces  are  unidirectional,  whereas  in  many  places  of  gray-matter, 
data  suggest  more  complicated  fiber  orientations.  The  partial  volume  effect  is 
apparent  on  the  left  image  where  the  corpus  callosum  neighbours  the  cortical  fibers. 

We  have  applied  the  same  procedure  to  the  image  of  an  excised  rat  spinal  cord 
acquired  at  14. IT.  The  results  can  be  seen  in  Figure  5-6.  The  rank  of  the  tensor 
model  employed  was  6.  Fiber  crossings  are  visible  in  many  areas  in  the  spinal  cord, 
particularly  in  the  ventral  nerve  roots  that  travel  among  white-matter  fiber  bundles 
in  a direction  perpendicular  to  them  and  causing  partial  volume  effects.  Also  note 
the  complicated  structure  in  gray-matter  where  most  fibers  are  oriented  in  the 
plane  of  the  image. 


In  the  previous  section  we  were  primarily  interested  in  the  fiber  orientations 
as  implied  by  a higher  rank  tensor.  We  have  done  this  by  generating  the  signal 
values  on  an  artificial  g-space  lattice,  and  then  taking  the  Fourier  transform 
numerically.  This  process  was  followed  by  an  isosurface  extraction  scheme.  All 
of  these  steps  take  a very  long  computational  time,  since  these  calculations  are 
repeated  for  each  voxel  of  the  image.  In  this  section  we  show  that  starting  from  the 
signal  attenuations  from  a IIARDI  experiment,  it  may  be  possible  to  calculate  the 
orientation  maps  without  the  need  to  fit  a high  rank  tensor  model  to  data. 

The  Fourier  transform  that  relates  the  signal  attenuation  to  the  water  dis- 
placement probability  (Eq.  2.21)  can  be  written  in  spherical  coordinates.  This  is  a 
consequence  of  the  pointwise  convergent  expansion  of  the  plane  wave  in  spherical 
coordinates  [87]  given  by 


5.4  A More  Direct  Approach  to  Orientation  Mapping 
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where  q = gg  and  AR  = Rr.  Inserting  this  expression  into  Eq.  2.51,  we  get 


Table  5 1:  Ai( g)  and  /?/( g)  functions  up  to  l = 6.  In  this  table,  (3  stands  for  /3(g). 


l 

Ms) 

Ms) 

0 

1 

0 

2 

-(i  + er2) 

3 

4 

(1  + 20/3-2  + 210/3-4) 

f (1  - 14/3-2) 

6 

- (1  + 42/3-2  + lf^/3-4  + 10395/3- 6) 

(1  - 36/?-2  + 396/3“4) 

where 

POO 

/;( g):=in  dq q2 ji(2irqRQ)  exp(— 47r2g2fD(g))  . (5.6) 

3o 

Note  that  the  function  Ps(Ro r)  is  not  the  isosurface  of  the  propagator  like  we 
attempted  to  construct  in  the  previous  section,  but  it  is  the  probability  of  finding 
the  particle,  initially  at  the  origin,  at  the  point  P0r,  i.e. , we  will  be  interested  in 
the  probability  values  on  a sphere  of  radius  R0.  The  integral  in  Eq.  5.6  can  be 
evaluated  and  it  is  given  by 
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where  iF±  is  the  confluent  hypergeometric  function.  These  functions  can  be  written 
as  the  sum  of  a Gaussian  and  an  error  function,  i.e., 
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The  kl/(g)  and  Bi( g)  functions  up  to  l = 6 are  given  in  Table  5-1. 

The  expression  in  Eq.  5.5  is  not  a convergent  expansion  of  P(i?Qr)  over  the 
index  l.  In  fact,  for  large  values  of  (3  the  higher  order  terms  can  be  larger  than 
the  lower  order  terms.  This  fact  can  be  seen  in  Figure  5-7.  Moreover,  using  the 
asymptotic  expansion  of  the  confluent  hypergeometric  functions  [88],  it  can  be  seen 
that  the  asymptotic  expansion  of  //(g)  functions  as  (3  — > oo  are  given  by 
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Figure  5-7:  Dependence  of  the  radial  integral  R on  (3.  The  curve  is  drawn  for  l 
values  from  0 to  6. 

Here,  the  first  term  is  subdominant  to  the  second  term  and  the  remainder  after  the 
truncation  at  an  order  l — L does  not  tend  to  0 as  /3  — > oo.  Therefore,  the  series 
in  Eq.  5.5  is  not  an  asymptotic  series  [89]. 

However,  we  speculate  that  the  angular  structure  of  the  probability  on  a 
sphere  of  radius  R0  may  be  retained  in  the  first  few  terms  of  the  expansion  if  a 
reasonable  choice  of  hence  (3  is  made.  In  Figure  5-  7,  we  see  that  for  large 
values  of  /?,  higher  order  terms  become  more  pronounced  relative  to  low  order  ones. 
Also,  in  the  proximity  of  the  origin,  the  higher  order  terms  have  a very  sharp  peak. 
Therefore,  a reasonable  choice  of  /3  will  be  between  these  two  extreme  behaviors. 
The  curves  in  Figure  5-  7 were  calculated  with  Dt  = .015mm2.  Evidently,  in  this 
case  a (3  value  of  1 to  3 may  yield  satisfactory  results.  We  have  tried  the  same 
procedure  for  different  values  of  diffusion  coefficients  observed  in  real  data,  and 
found  that  this  choice  of  acceptable  range  for  (3  has  a weak  dependence  on  the 
diffusivities. 
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In  order  to  estimate  the  probability  on  the  surface  of  a sphere  of  radius  Rq,  we 


go  back  to  Eqs.  5.5  and  5.6,  and  expand  //(g)  in  Laplace  series,  i.e., 
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Comparing  the  integration  over  g in  Eq.  5.5  with  the  expression  in  Eq.  5.12,  it  is 
easy  to  see  that 


which  is  just  a Laplace  series  expansion  of  Ps{Ro r).  Note  that  coefficients  of  this 
Laplace  series  for  some  l value  come  from  the  Z-th  order  Laplace  series  coefficients 


In  summary,  the  estimation  of  the  probability  of  finding  the  particle  at  the 
point  Ror  away  from  the  origin  involves  the  following  steps: 

• Given  HARDI  data  including  Nq  images  acquired  with  gradient  directions 
along  different  orientations,  calculate  the  diffusivity  D( g)  along  each  direc- 
tion. 

• Calculate  //(g)  using  Eq.  5.7  or  5.8  and  Table  5-1. 

• Calculate  ct//m,  the  l- th  order  spherical  harmonic  transform  of  //(g)  for  each  l. 

• Evaluate  Eq.  5.13,  which  is  just 
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5.4.1  Results 


We  have  applied  the  scheme  described  above  to  the  simulations  that  we  have 
performed  in  the  previous  section.  Figure  5-8  shows  the  effect  of  selection  for  /?  on 
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Figure  5-8:  Probability  maps  estimated  on  a sphere  of  radius  Rq,  which  is  deter- 
mined by  the  choice  of  (3  through  Eq.  5.9.  These  surfaces  were  constructed  by 
setting  (3  to  a value  of  from  0.75  to  4 in  equal  steps  (from  left  to  right).  Top  row 
shows  these  surfaces  when  there  is  only  one  orientation,  where  the  bottom  row 
shows  the  same  surfaces  from  a voxel  with  two  distinct  orientations. 


the  constructed  probability  surfaces.  Increasing  (3  results  in  sharpening  similar  to 
increasing  P0  value  when  the  isosurfaces  of  the  probability  function  was  calculated 
in  generalized  DTI  based  orientation  mapping.  Small  /3  forces  R0  to  be  smaller 
than  the  characteristic  length  y/Dt,  in  which  case  the  distribution  of  probability  on 
the  surface  becomes  more  uniform. 

We  have  also  calculated  the  probability  surfaces  for  the  simulated  image 
originally  introduced  in  Figure  4 1.  The  surfaces  are  very  consistent  with  the 
underlying  known  fibrous  structure.  It  is  also  consistent  with  the  generalized  DTI 
based  approach  we  presented  in  the  last  section. 

The  probability  maps  calculated  from  the  imaging  data  are  also  consistent 
with  the  previous  approach.  We  provide  in  Figure  5-10,  the  orientation  maps 
calculated  from  the  slices  that  were  presented  in  the  previous  chapter  in  Figure 
4-6.  In  that  section  our  goal  was  to  understand  whether  the  differences  observed 
in  the  anisotropy  and  variance  values  were  due  to  fiber  crossings.  Looking  at  those 
regions  where  we  had  predicted  complicated  structure  based  upon  the  anisotropy 
and  variance  changes,  we  see  in  Figure  5-4  that  there  are  indeed  orientational 
heterogeneity  in  those  voxels. 
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Figure  5 9:  The  probability  maps  from  the  simulated  image  of  Figure  4 1 calcu- 
lated using  the  expansion  of  the  probability  on  the  surface  of  a sphere. 

5.5  Discussion 

Prior  to  DTI,  diffusion  weighted  imaging  (DWI)  had  been  implemented  and 
shown  to  be  very  sensitive  to  many  pathologies,  where  the  most  striking  result 
was  in  the  diagnosis  of  cerebral  ischemia  [25].  However  in  DWI,  the  oriented 
structure  of  the  neural  tissue  was  a problem  since  the  choice  of  diffusion  gradient 
direction  could  change  the  outcome  of  the  measurement.  DTI  exploited  this  very 
phenomenon  by  proposing  the  rank-2  tensor  as  a model  that  could  quantify  the 
orientational  dependence  of  the  measurement,  and  based  on  this,  it  produced 
meaningful  results  in  quantifying  the  level  of  anisotropy  [57]  and  determining  the 
fiber  directions  in  major  white-matter  regions  [47,  90].  The  quantity  that  could  be 
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Figure  5-10:  The  probability  maps  calculated  in  the  selected  slices  of  the  rat  brain. 
These  slices  are  the  same  with  those  given  in  Figure  4-6.  The  region  from  which 
these  surfaces  are  depicted  in  blue  in  the  GA  maps  provided  in  the  lower  left  corner 
of  each  image. 

quantified  by  DWI  was  the  apparent  diffusion  coefficient,  which  is  a rank-0  tensor. 
In  this  context,  the  transition  from  quantitative  DWI  to  DTI  is  an  increase  in  the 
rank  of  the  tensor  being  evaluated.  We  have  presented  a novel  approach  to  extend 
this  to  higher  rank  tensors  to  overcome  some  of  the  difficulties  experienced  with 
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traditional  (rank-2)  DTI.  We  showed  that  it  is  possible  to  map  the  orientations 
more  accurately  using  these  higher  rank  tensors. 

Traditional  DTI  assumes  an  oriented  Gaussian  displacement  profile,  which  is 
true  under  the  assumption  that  the  signal  attenuation  along  an  arbitrary  direction 
will  be  monoexponential.  This  is  known  to  be  not  the  case  in  the  presence  of 
restrictions.  However,  the  success  of  DTI  in  determining  the  fiber  orientations  in 
major  white-matter  tracts  suggests  that  this  assumption  is  a good  one  in  the  sense 
that  assuming  monoexponentiality  doesn’t  affect  the  fiber  orientations  significantly. 
When  estimating  the  displacement  probabilities,  our  approach  uses  the  same 
assumption  but  for  more  complicated  diffusivity  profiles  than  those  that  can  be 
expressed  by  traditional  (rank-2)  DTI.  The  simulation  results  further  verified  that 
monoexponentiality  did  not  have  an  affect  on  the  fiber  orientations  observed  at 
least  for  relatively  simple  cases.  Note  that  since  the  diffusivity  and  probability 
profiles  have  the  same  orientational  behavior  in  one-fiber  case,  it  was  sufficient 
in  the  case  of  rank-2  DTI  to  diagonalize  the  tensor  and  calculate  the  principal 
eigenvector  to  find  the  direction  of  the  fiber.  Since  diffusivity  profiles  in  general  do 
not  give  the  fiber  directions,  when  a higher  rank  tensor  model  is  used,  one  has  to 
calculate  the  displacement  profiles  to  resolve  more  complicated  structures. 

In  the  latter  part  of  this  chapter,  we  presented  a more  direct  approach  to 
probability  mapping  that  bypassed  the  calculation  of  a tensor  from  the  signal 
attenuation  profiles.  This  scheme  employed  the  same  monoexponentiality  assump- 
tion that  was  made  in  traditional  as  well  as  generalized  DTI.  The  Laplace  series 
coefficients  of  the  radial  integral  involving  the  signal  attenuations  were  related  to 
the  Laplace  series  coefficients  of  the  probability  on  the  sphere.  We  kept  only  a few 
terms  in  this  Laplace  series  although  we  were  unable  to  estimate  the  error  involved 
in  this  process.  Nevertheless  the  performance  of  the  algorithm  was  highly  satisfac- 
tory both  in  simulations  and  real  (noisy)  data.  This  is  most  likely  because  although 
the  calculated  values  may  be  wrong,  the  angular  structure  of  these  values  might 
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be  quite  accurateeven  when  only  a few  terms  are  included  in  the  expansion.  A 
significant  gain  that  was  achieved  by  using  this  procedure  rather  than  the  previous 
one  employing  higher  rank  tensors  was  the  computational  speed.  We  have  found 
that,  although  we  have  not  spent  any  extra  effort  in  making  the  algorithm  efficient, 
this  scheme  was  about  40-50  times  faster  than  the  previous  approach. 

A well-known  method  commonly  used  to  observe  porous  structures  is  the 
g-space  imaging  method  [54],  which  was  also  proposed  for  use  in  finding  the 
orientations  in  tissue  [91].  In  this  approach,  an  image  is  acquired  at  every  point 
of  the  Cartesian  g-space  grid.  Even  for  a very  sparse  sampling  such  as  on  an  8 x 
8x8  grid,  a total  of  512  acquisitions  is  necessary.  Note  that  to  avoid  truncation 
effects  one  has  to  go  into  considerable  depth  in  g-space,  in  which  case  the  acquired 
images  will  be  very  low  in  signal-to-noise  ratio  because  of  the  very  large  signal 
attenuation.  Although  g-space  imaging  is  model  independent  in  theory,  feasibility 
of  covering  only  a relatively  small  grid  limits  the  extent  to  which  the  structure 
in  tissue  can  be  resolved.  The  approaches  we  have  presented  on  the  other  hand 
require  fewer  number  of  acquisitions  where  all  of  these  are  along  different  directions 
in  g-space.  Therefore,  the  emphasis  is  solely  on  extracting  the  orientational 
information  from  the  sample.  Collecting  fewer  samples  (in  total)  at  lower  6-values 
makes  the  technique  less  demanding  in  acquisition  time  as  well  as  in  terms  of  the 
gradient  coils  to  be  used.  In  our  simulations,  we  found  it  satisfactory  to  use  data 
with  diffusion  gradients  applied  along  46  directions.  However,  larger  numbers  are 
always  beneficial  especially  in  the  presence  of  noise.  As  a result,  for  the  excised  rat 
brain  data  that  we  have  shown,  although  we  had  81  measurements,  which  could 
accomodate  a rank-10  model  in  our  first  approach,  we  chose  to  use  the  rank-8 
model  so  that  the  sensitivity  of  the  results  to  noise  is  reduced.  Similarly,  truncation 
of  the  series  at  l = 6 in  the  second  scheme  served  the  same  purpose.  When  this  is 
done,  voxels  with  even  higher  complexity  are  likely  to  be  seen  as  isotropic. 


CHAPTER  6 

DISCUSSION  AND  CONCLUSIONS 

In  this  dissertation  we  presented  new  techniques  that  enable  one  to  construct 
more  accurate  anisotropy  and  orientation  maps  compared  to  those  achieved  by 
the  commonly  employed  DT-MRI  approach.  The  main  motivation  was  to  be  able 
to  enhance  accuracy  without  imposing  requirements  that  are  not  easy  to  attain 
in  clinical  conditions.  Although  we  have  not  used  data  collected  in  a clinical 
setting,  the  diffusion  weighting  factors  and  the  number  of  images  required  is  not 
too  demanding  for  it  to  be  extended  to  the  clinical  environment. 

Our  approach  in  the  generalization  of  DT-MRI  was  to  extend  the  transition 
from  diffusion  weighted  to  diffusion  tensor  imaging  by  incorporating  tensors  of 
rank  higher  than  2 into  the  formalism  relating  the  signal  attenuations  to  the 
diffusional  properties  of  the  sample.  This  extension  reduced  the  loss  of  information 
encountered  in  the  fitting  of  the  rank-2  tensor  while  using  slightly  overdetermined 
data  had  the  effect  of  reducing  the  sensitivity  of  the  approach  to  noise. 

Although  fiber  orientation/tract  mapping  is  the  primary  goal  of  the  studies  on 
the  measurement  and  modeling  of  anisotropic  diffusion  in  fibrous  tissue,  most  of  the 
clinical  studies  quantify  the  multiorientational  images  using  rotationally  invariant 
scalar  measures  derived  from  these  images.  Because  of  the  widespread  utilization 
of  scalar  indices  in  clinically  oriented  studies,  we  discussed  the  construction  of 
scalar  measures  such  as  mean  diffusivity  and  anisotropy  and  suggested  several 
other  measures  that  can  be  useful  in  the  upcoming  studies.  We  have  shown  that 
the  difficulties  associated  with  DT-MRI  creates  problems  in  the  anisotropy  maps 
which  may  have  significant  influence  on  these  clinical  studies.  The  scalar  indices  we 
have  presented  are  not  designed  for  a particular  rank  tensor  model,  but  it  can  be 
utilized  with  tensors  of  arbitrary  rank  and  even  functions  defined  on  the  surface  of 


98 


99 

the  sphere.  This  makes  it  possible  to  make  meaningful  comparisons  of  the  values  of 
these  indices  calculated  through  different  methods. 

The  employment  of  higher  rank  tensors  was  also  useful  in  creating  more 
accurate  orientation  maps.  Although  the  calculation  of  these  orientation  maps 
was  a computationally  difficult  task,  it  demonstrated  that  the  assumption  made 
in  these  calculations,  namely  the  monoexponential  decay  along  each  radial  line 
of  space  defined  by  the  diffusion  sensitizing  gradients,  was  acceptable  for  fiber 
orientation  mapping.  Using  the  same  assumption,  we  attempted  to  construct  the 
orientation  maps  in  a more  efficient  way  which  yielded  orientation  maps  of  similar 
quality. 

Although  we  constructed  fiber  orientation  maps  based  on  the  techniques  we 
have  developed,  we  did  not  produce  fiber-tract  maps  like  we  did  in  the  case  of  DT- 
MRI  because  when  there  are  several  voxels  with  multiple  orientations  in  a region, 
the  calculation  of  the  tracts  become  severely  ill-posed  in  contrast  with  DT-MRI 
that  is  deterministic  since  it  yields  only  one  orientation  at  each  point  in  space. 
Generating  fiber-tracts  from  images  with  multiorientational  voxels  remains  an  open 
problem. 

Throughout  our  studies  we  found  it  useful  to  test  the  techniques  we  developed 
on  simulated  voxels  and  synthetic  images.  These  simulations  served  as  a necessary 
condition  for  the  applicability  of  the  techniques  to  real  data.  The  knowledge  gained 
through  simulations  also  helped  us  interpret  the  results  from  real  data  correctly. 

Majority  of  the  studies  on  diffusion  MRI  of  neural  tissue  including  this 
dissertation  have  used  the  diffusion  gradient  strength  and  orientation  as  the  control 
variable  in  diffusion  experiments.  However,  the  diffusional  process  is  dependent 
upon  other  factors  that  can  also  be  controlled  in  MRI  experiments.  These  include 
diffusion  time  and  diffusion  gradient  pulse  length.  However,  incorporation  of  these 
parameters  into  the  formalism  of  restricted  diffusion  is  a challenging  and  open 
problem. 
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Diffusion  imaging  of  the  neural  tissue  holds  great  promise  in  the  inference 
of  other  information  about  tissue  microstructure.  For  example,  by  observing 
and  modelling  restricted  diffusion,  it  may  be  possible  in  the  future  to  quantify 
parameters  such  as  surface  to  volume  ratio,  mean  curvature  of  the  boundaries, 
and  permeability  maps.  These  parameters  are  very  significant  especially  because 
they  are  altered  by  the  pathological  perturbations  associated  with  many  cerebral 
diseases. 


APPENDIX 

ACQUISITION  PARAMETERS 

The  images  we  have  shown  throughout  this  dissertation  are  from  data  acquired 
at  the  Advanced  Magnetic  Resonance  Imaging  and  Spectroscopy  Facility  (a  section 
of  the  National  High  Magnetic  Field  Laboratory)  at  McKnight  Brain  Institute 
of  University  of  Florida.  In  this  appendix  we  provide  the  parameters  of  the 
acquisitions  we  have  performed  on  several  tissue  samples.  The  acquisitions  are 
presented  in  order  of  appeareance  in  the  text.  All  experiments  were  performed 
with  the  approval  of  the  University  of  Florida  Institutional  Animal  Care  and  Use 
Committee.  All  tissue  was  fixed  in  4%  paraformaldehyde  and  imaged  in  phosphate 
buffer  saline. 

Data  from  excised  rat  brain  presented  in  Figures  2-3,  2-5,  2-12  and 
2-13.  An  excised  rat  brain  was  imaged  at  17. 6T  (750 MHz  proton  frequency) 
using  a Bruker  Avance  imaging  spectrometer.  A diffusion  weighted  spin  echo  pulse 
sequence  was  used  with  diffusion  gradient  strengths  of  335,  475  and  585 mT/m 
producing  b- values  of  approximately  750,  1500,  2250 s/mm2.  Diffusion  gradients 
were  applied  along  seven  directions:  ex,  ey,  e2,  -~^(ex  + ey),  -^(ex  + ez),  -T(ey  + ez), 
^(ex  + ey  + e2).  Measurement  in  each  direction  and  b- value  was  repeated  6 times. 
Imaging  parameters  were:  TR  = 3000ms,  TE  = 27.7ms,  A = 17.8ms,  8 — 2.4ms, 
resolution=  117  x 117  x 250/rm3,  held  of  view=  15  x 15  x 19.5mm3  and  the  matrix 
size  was  128  x 128  x 78. 

Data  from  excised  rat  hippocampus  presented  in  Figures  2-6  and  2— 

7.  An  excised  rat  hippocampus  was  imaged  at  14.  IT  (600 MHz  proton  frequency) 
using  a Bruker  Avance  imaging  spectrometer.  A diffusion  weighted  stimulated 
echo  pulse  sequence  was  used  with  a total  of  17  diffusion  gradient  strengths  from 
0 to  2936 mT/m  in  steps  of  183.5 mT/m  all  of  which  were  applied  along  the  main 
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magnetic  field.  The  corresponding  b- values  varied  approximately  from  0 to, 
50000s/mm2.  Each  measurement  was  repeated  12  times.  Imaging  parameters  were: 
TR  = 1000ms,  TE  = 12.6ms,  A = 20ms,  5 = 2ms,  resolution=  78  x 78  x 750/rm3, 
field  of  view=  5 x 5 x 1.5mm3  and  the  matrix  size  was  64  x 64  x 3. 

Data  from  excised  rat  spinal  cord  presented  in  Figures  2-9,  2-10, 
2—11,  2-12  and  2-13.  An  excised  rat  spinal  cord  was  imaged  at  17. 6T  (750 MHz 
proton  frequency)  using  a Bruker  Avance  imaging  spectrometer.  A diffusion 
weighted  spin  echo  pulse  sequence  was  used  with  diffusion  gradient  strengths 
of  0,  180,  360  and  540 mT/m  producing  b- values  of  approximately  18,  245,  925, 
2060s/mm2.  Diffusion  gradients  were  applied  along  seven  directions:  ex,  ey,  e2, 
^j(ex+ey),  ^=(ex+e2),  ^(ey+e2),  ^(ex+ey+e2).  Measurement  in  each  direction 
and  6-value  was  repeated  16  times.  Imaging  parameters  were:  TR  = 1000ms, 

TE  = 27.1ms,  A = 17.8ms,  6 = 2.4ms,  resolution=  40  x 40  x 250/im3,  field  of 
view=  5.1  x 5.1  x 5mm3  and  the  matrix  size  was  128  x 128  x 20. 

Data  from  excised  normal  and  injured  rat  spinal  cords  presented  in 
Figure2-14.  Both  spinal  cords  were  imaged  at  17.6T  (750 MHz  proton  frequency) 
using  a Bruker  Avance  imaging  spectrometer.  A diffusion  weighted  spin  echo  pulse 
sequence  was  used  with  diffusion  gradient  strengths  of  112,  264  and  375 mT/m 
producing  b- values  of  approximately  100,  500,  1000 s/mm2.  Diffusion  gradients 
were  applied  along  seven  directions:  ex,  ey,  e2,  ^(ex  + ey),  ^(ex  + ez),  ^(ey  + e2), 
^(ex  + ey  + e2).  Measurement  in  each  direction  and  6-value  was  repeated  13  times. 
Imaging  parameters  were:  TR  = 1570ms,  TE  — 28.8ms,  A - 17.8ms,  5 = 2.4ms, 
resolution=  39  x 39  x 250/rm3,  field  of  view=  5 x 5 x 10mm3  and  the  matrix  size 
was  128  x 128  x 40. 

Data  from  excised  rat  brain  presented  in  Figure  2-15.  An  excised  rat 
brain  was  imaged  at  17. 6T  (750 MHz  proton  frequency)  using  a Bruker  Avance 
imaging  spectrometer.  A diffusion  weighted  spin  echo  pulse  sequence  was  used 
with  diffusion  gradient  strengths  of  119,  267  and  377 mT/m  producing  b-values 
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of  approximately  100,  500,  1000 s/mm2.  Diffusion  gradients  were  applied  along 
seven  directions:  ex,  e„,  ez,  -^(ex  + ey),  ^(ex  + e2),  ^(ey  + e2),  ^(ex  + ey  + 
e2).  Measurement  in  each  direction  and  b- value  was  repeated  9 times.  Imaging 
parameters  were:  TR  = 3060ms,  TE  = 28.8ms,  A = 17.8ms,  5 = 2.4ms, 
resolution=  117  x 117  x 270/rm3,  field  of  view=  15  x 15  x 21mm3  and  the  matrix 
size  was  128  x 128  x 78. 

Data  from  excised  rat  brain  presented  in  Chapter  3.  An  excised 
rat  brain  was  imaged  at  11.  IT  (470 MHz  proton  frequency)  using  a Bruker 
A vance  imaging  system.  A diffusion  weighted  spin  echo  pulse  sequence  was  used. 
Diffusion  weighted  images  were  acquired  along  81  directions  specified  by  the 
third  order  tessellations  of  an  icosahedron  on  the  hemisphere  with  a gradient 
strength  of  150 mT/m,  producing  a b- value  of  1050mm2/s  along  with  a single  image 
acquired  at  b ss  0.  Number  of  repetitions  was  4.  The  imaging  parameters  were 
TE  = 32.9ms,  TR  = 2350ms,  A = 20.0  ms,  5 = 6.0  ms,  resolution=  187.5  x 187.5  x 
1000/xm3,  field  of  view=  2.4  x 2.4  x 24mm3  and  the  matrix  size  was  128  x 128  x 24. 

Data  from  excised  rat  brain  presented  in  Chapters  4 and  5.  An  excised 
rat  brain  was  imaged  at  17.6T  (750 MHz  proton  frequency)  using  a Bruker  Avance 
imaging  system.  A diffusion  weighted  spin  echo  pulse  sequence  was  used.  Diffusion 
weighted  images  were  acquired  along  81  directions  specified  by  the  third  order 
tessellations  of  an  icosahedron  on  the  hemisphere  with  a gradient  strength  of 
500 mT/m,  producing  a b- value  of  1500mm2/s  along  with  a single  image  acquired 
at  b « 0.  Number  of  repetitions  was  6.  The  imaging  parameters  were  TE  = 28ms, 
TR  — 2500ms,  A = 17.8ms,  8 = 2.2ms,  resolution=  150  x 150  x 300/xm3,  field  of 
view=  15  x 15  x 18mm3  and  the  matrix  size  was  100  x 100  x 60. 

Data  from  excised  rat  spinal  cord  presented  in  Chapter  5.  An  excised 
rat  brain  was  imaged  at  14. IT  (600 MHz  proton  frequency)  using  a Bruker  Avance 
imaging  system.  A diffusion  weighted  spin  echo  pulse  sequence  was  used.  Diffusion 
weighted  images  were  acquired  along  46  directions  specified  by  the  second  order 
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tessellations  of  an  icosahedron  on  the  hemisphere  with  a gradient  strength  of 
733 mT/m,  producing  a 6- value  of  1500mm2 / s along  with  a single  image  acquired 
at  b 0.  Number  of  repetitions  was  7.  The  imaging  parameters  were  TE  = 25 ms, 
TR  = 1170ms,  A = 17.5ms,  6 = 1.5ms,  resolution=  59.7  x 59.7  x 300/xm3,  field  of 
view=  4.3  x 4.3  x 12mm3  and  the  matrix  size  was  72  x 72  x 40. 
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